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This  dissertation  is  a  theoretical  investigation  of  col 
lisional  dynamics  of  two  particle  interactions  in  a  gravita¬ 
tional  field.  This  research  is  unique  in  that  it  is  the 
first  attempt  at  modeling  the  hydrodynamic  interactions  be¬ 
tween  a  nonspherical  particle  and  a  spherical  particle  under 
going  gravitational  collisions  in  an  LMFBR  environment. 
First,  basic  definitions  and  expressions  are  developed  for 
nonspherical  particles  and  related  to  spherical  particles  by 
means  of  shape  factors.  Using  volume  equivalent  diameter  as 
the  defining  length  in  the  gravitational  collision  kernel, 
the  aerodynamic  shape  factor,  k,  the  density  correction  fac¬ 
tor,  a,  and  the  gravitational  collision  shape  factor,  B,  are 
used  to  correct  the  collision  kernel  for  the  case  of  colli¬ 
sions  between  aerosol  agglomerates.  The  Navier-Stokes  equa¬ 
tion  in  oblate  spheroidal  coordinates  is  solved  to  model  a 
nonspherical  particle  and  then  the  dynamic  equations  for  two 
particle  motions  are  developed.  A  computer  program  NGCEFF 
is  constructed,  the  Navier-Stokes  equation  is  solved  by  the 
finite  difference  method,  and  the  dynamical  equations  are 
solved  by  Gear's  method.  Results  are  obtained  for  several 


cases  and  are  compared  with  previous  work  of  atmospheric 
sciences  and  LMFBR  studies  for  the  spherical  cases.  It  is 
concluded  that  the  aerosol  gravitational  collision  shape 
factor  can  be  determined  by  further  theoretical  work  based 
on  the  concepts  and  methods  developed  in  this  dissertation. 
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CHAPTER  I 


INTRODUCTION 


1 . 1  LMFBR  Safety  Requirements 

Safety  analyses  for  liquid  metal  fast  breeder  reactors 
(LMFBR)  are  based  on  predictions  of  radioactivity  released 
to  the  environment  under  a  postulated  serious  accident  con¬ 
dition  such  as  a  Hypothetical  Core  Disruptive  Accident 
(HCDA)  which  presupposes  fuel  meltdown  and  vaporization,  and 
the  release  of  coolant,  fuel  and  structural  material  into 
the  primary  containment.  In  view  of  the  light  water  reactor 
licensing  requirements,  it  is  clear  that  detailed  analyses 
of  the  consequences  of  an  HCDA  will  be  needed.  A  critical 
step  in  these  analyses  for  LMFBR  is  the  prediction  of  mass 
concentration  for  radioactive  material  airborne  in  the  con¬ 
tainment  building  at  various  times  after  the  postulated 
accident.  If  this  can  be  accurately  calculated,  the  charac¬ 
teristics  of  radioactive  aerosols  outside  the  containment 
vessel  can  be  predicted  as  a  function  of  time.  Armed  with 
this  information,  general  population  exposure  levels  down¬ 
wind  of  the  reactor  site  can  be  estimated.  Consequently, 
the  development  of  accurate  and  physically  realistic  aerosol 
behavior  models  has  been  the  subject  of  extensive  studies. 
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To  help  visualize  the  HCDA  environment.  Figure  1.1.1 
is  a  schematic  of  the  proposed  Clinch  River  Breeder  Reactor 
(CRBR)  containment  building.  The  accident  initiating  events 
presently  have  not  been  defined.  Only  the  damage  mechanisms 
have  been  studied*.  Starting  with  an  energetic  power  excur¬ 
sion,  the  active  portion  of  the  core  is  heated  to  a  uniform 
temperature  of  4800°K.  The  core  material  expands  upward  and 
penetrates  the  sodium  pool  as  a  hemispherical  bubble  com¬ 
posed  of  solid  fragments,  aerosols  droplets,  vapor,  etc. 
(Figure  1.1.2).  The  pressure  wave  created  by  the  bubble 
damages  the  reactor  vessel,  the  closure  head,  and  the  guard 
vessel.  In  most  postulated  scenarios,  the  sodium  pool  burns 
and  creates  a  log-normal  aerosol  size  distribution  with  a 

3 

maximum  concentration  not  exceeding  200  g/m  .  The  aerosol 
inside  the  containment  building  is  considered  well-mixed 
with  no  spatial  inhomogeneities  for  the  aerosol. 

Thus,  following  an  HCDA,  an  aerosol  mixture  of  sodium 
compounds,  fuel  and  core  structural  materials  will  begin  to 
plate  out  on  walls,  floors,  equipment  and  other  components 
because  of  various  mechanisms.  One  key  mechanism  is  gravi¬ 
tational  coagulation  between  aerosol  particles  undergoing 
simple  gravitational  settling.  This  affects  the  particle 
size  distribution  and  hence  the  release  of  radioactivity  to 
the  outside  environment. 
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Various  models  have  been  developed  to  describe  aerosol 
behavior  under  postulated  accident  scenarios.  Examples  of 
such  models  are  the  HAARM-3  (Heterogeneous  Aerosol  Agglomer¬ 
ation  Revised  Model  -  3)  developed  at  Battelle  Columbus 
Laboratories^,  the  PARDISEKO- 1 1 1  code  at  the  Karlsruhe 

Nuclear  Research  Center^  and  the  TRAP  code  by  Brookhaven 

4 

National  Laboratory  .  These  codes  attempt  to  solve  the 
aerosol  rate  equation,  an  integro-differential  equation 
describing  the  rate  of  change  of  particle  concentration  due 
to  the  various  agglomeration  and  deposition  mechanism. 

1 . 2  Dissertation  Organization 

This  research  investigates  the  hydrodynamic  interaction 
between  two  LMFBR  aerosols  falling  in  a  gravitational  field. 
The  probability  for  particle-particle  collision,  due  only  to 
the  hydrodynamic  and  gravitational  forces,  can  be  described 
by  defining  the  gravitational  collision  efficiency  of  inter¬ 
acting  particles.  Determined  experimentally  or  theoretically 
the  collision  efficiency  is  used  to  estimate  the  growth  of 
aerosol  particles. 

The  dissertation  investigates  several  ways  to  improve 
current  theoretical  methods  for  determining  the  gravitational 
collision  efficiency.  Chapter  II  is  a  review  of  LMFBR 
aerosol  characteristics,  such  as  size  distribution,  chemical 
composition,  shape  and  density.  The  aerosol  rate  equation 
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is  presented  and  the  gravitational  collision  cross  section 
and  efficiency  are  defined. 

Chapter  III  contains  a  review  of  previous  work  done  to 
estimate  the  gravitational  collision  efficiency.  Determining 
the  gravitational  collision  kernel  is  a  fundamental  problem 
of  aerosol  science  which  has  been  extensively  investigated 
by  meteorologists  because  of  its  importance  in  the  formation 
and  growth  of  raindrops  and  hailstones. 

Chapter  IV  derives  the  equations  needed  to  calculate  a 
generalized  gravitational  collision  kernel.  A  review  of  the 
HAARM- 3  code  is  presented  using  the  expressions  derived  from 
this  chapter. 

Chapter  V  presents  the  numerical  methods  used  to  solve 
the  equations  and  the  problems  associated  with  various  nu¬ 
merical  techniques.  Code  development  is  examined  and  options 
available  are  explained. 

Chapter  VI  contains  the  results  for  the  cases  studied 
and  their  application  to  current  aerosol  codes  like  HAARM-3. 

Chapter  VII  reviews  the  research  project  and  makes 
recommendations  for  future  work. 

Studies  associated  with  the  dissertation  topic  are 
presented  in  the  appendices.  It  includes  kinetic  corrections 
to  the  superposition  method  for  calculating  the  gravitational 
collision  kernel  for  spherical  particles,  a  listing  of 
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FORTRAN  function  routines  for  calculating  gravitational 
collision  efficiencies,  and  a  discussion  of  the  computer 
codes  used  for  the  dissertation. 


CHAPTER  II 


LMFBR  AEROSOLS  AND  GOVERNING  EQUATIONS 

2 . 1  Morphology  and  Aerodynamics  of  LMFBR  Aerosols 

Aerosols  can  be  produced  by  fragmentation  during  the 
postulated  energetic  core  disruptive  accident  and  by  conden¬ 
sation  of  vapors  during  and  following  the  accident.  The 
resulting  aerosols  can  consist  of  a  number  of  different 
materials,  such  as  sodium,  plutonium,  uranium  and  their 
oxides,  and  reactor  structural  material.  Due  to  the  large 
mass  of  sodium  coolant  present,  it  is  reasonable  to  assume 
that  the  bulk  of  LMFBR  aerosols  is  composed  of  sodium  com¬ 
pounds.  Because  of  this  assumption,  the  majority  of  investi¬ 
gators  have  reported  the  characteristics  and  behavior  of 
sodium  and  its  compounds,  although  recently  more  research 
into  the  morphology  of  plutonium  and  uranium  oxides  has  been 
reported. 

2.1.1  Primary  Particles 

Condensation  is  believed  to  be  more  important  to  the 
production  of  LMFBR  aerosols  than  fragmentation^.  Particles 
formed  by  nucleation  are  of  the  order  of  0.001  ym  diamter. 
Once  formed,  they  grow  rapidly  by  condensation.  Although 
primary  particles  actually  define  the  starting  point  for  the 
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agglomeration  process,  they  coagulate  so  rapidly  that  it 
is  impractical  to  describe  the  aerosol  source  in  terms  of 
the  characteristics  of  primaries.  Therefore,  the  aerosol 
term  is  normally  characterized  by  specifying  the  properties 
of  the  LMFBR  aggregates. 

2.1.2  LMFBR  Aggregates 


It  is  generally  accepted  that  the  agglomerates  of  sodium 
oxide,  plutonium  and  uranium  oxides  are  branched  chain-like 
structures  composed  of  spherical,  small  in  diameter  (0.01  to 
0.45  urn)  primary  particles  that  are  distributed  log-nor- 

ft  -  ft 

mally  .  Because  of  the  complex  nature  of  the  structures, 
measurements  of  these  agglomerates  are  typically  concerned 
with  aerodynamic  shape  factor,  k,  and  density  correction 
factor,  a,  given  mathematically  by  the  formulas: 


(2. 1.2.1) 


a  =  —  (2. 1.2. 2) 

pm 

where:  V$  ~  Stokes  settling  velocity 

Va  ~  Actual  or  measured  velocity 
Pa  -  Density  of  the  agglomerate 
Pm  -  Density  of  the  material. 


10 

These  quantities  are  further  defined  in  Chapter  IV. 

The  important  point  is  that  solid  aerosol  particles  are 
rarely  seen  to  be  spherical,  at  least  so  far  as  those  par¬ 
ticles  are  concerned  which  can  be  seen  in  an  optical  micro¬ 
scope  or  electron  microscope.  Figures  2.1. 2.1  and  2. 1.2. 2 
typify  the  general  characteristics  of  chain  agglomerates  and 
indicate  the  problem  associated  with  any  mathematical 
description  of  the  aggregates.  Aggregates  have  been  classi¬ 
fied  into  two  groups,  one  appearing  to  be  almost  clusters 
of  particles  and  the  second  type  being  long  chains.  The 
theoretical  problem  is  to  simulate  these  agglomerates  with 
simple  geometric  solids  which  characterize  the  general 
behavior  of  the  aerosol. 

Jordan  and  Geiseke6  measured  the  aerodynamic  shape 
factor  for  agglomerates  of  UC^  aerosol  particles  and  esti¬ 
mated  the  approximate  range  as  3  I  <  £  16.  Using  the  work 
9 

of  Kops  et  al .  ,  Jordan  and  Geiseke  tried  to  correlate  k 
with  mass  equivalent  radius  of  the  agglomerate  but  were 
unsuccessful.  Their  conclusion  was  that  U02  particles 
observed  in  the  Millikan  thermal  cell  did  not  fall  into  the 
category  of  either  fluffy  ball  with  constant  packing  density 
or  straight  chain.  This  conclusion  is  difficult  to  support 
because  the  statistics  of  the  Millikan  analysis  is  poor  com¬ 
pared  with  that  of  the  Stober  centrifuge  method,  which  was 
the  technique  used  by  Kops.  Jordan  and  Gieseke  measured  the 
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characteristics  of  individual  particles  in  the  Millikan  cell 
whereas  Kops  measured  the  average  properties  of  agglomerates 
with  the  Stober  centrifuge.  Kops'  work  on  the  aerodynamic 
diameter  of  branched  chain-like  agglomerates  was  verified  by 
Wegrzyn  and  Shaw^  who  used  a  sodium  burn  chamber  to  produce 
the  metal  oxide  instead  of  the  exploding  wire  technique  of 
Kops.  Their  results  seem  consistent  with  the  observed  fact 
that  the  metal  oxide  vapors  condense  to  form  aggregates. 

O 

The  most  consistent  work  is  that  of  Wegrzyn  and  Shaw  , 
Kops  et  al .  ,  and  Van  de  Vate  et  al .  .  These  studies  indi¬ 
cate  that  the  shape  of  the  metal  oxide  aerosol  is  a  branched, 
chain-like  agglomerate  made  up  of  spherical  primary  particles. 
Spiral  centrifuge  analysis  yielded  aerodynamics  shape  factors, 
k,  between  0.8  and  4.  Van  de  Vate  observed  that  increased 
humidity  for  sodium  oxide  particles  would  lead  to  compaction 
of  the  fluffy  aggregates  resulting  in  an  increase  aerodynamic 
diameter  and  decrease  shape  factor.  None  of  the  researchers 
reported  a  shape  factor  greater  than  4  and  thus  the  work  of 

Jordan  and  Gieseke  cannot  be  verified. 

12  18 
Pertmer  ,  and  Pertmer  and  Loyalka  showed  that  the 

gravitational  collision  efficiency  is  sensitive  to  the  par¬ 
ticle  density.  Pertmer  reported  Figure  2. 1.2. 3  showing  the 
relationship  between  agglomerate  radius  and  density.  This 
result  is  from  Battelle  Columbus  Laboratories  and  has  not 

O 

been  verified  elsewhere.  Wegrzyn  and  Shaw  did  similar 


sodium  pool  fire  studies  and  reported  an  average  agglomerate 
density  (apm)  of  1.32  g/cm3  for  chain-like  particles  with 
aerodynamic  diameters  between  0.6  to  0.7  um.  If  one  assumes 
a  primary  particle  density  of  2.27  g/cm  (Na^O)  then  the 
density  factor  (a)  is  0.58.  There  is  evidence^  that  as  the 
agglomerate  grows,  the  density  factor  decreases  as  one  would 
predict  but  a  value  as  low  as  0.13  (suggested  by  Figure 
2. 1.2. 3)  has  not  been  collaborated.  Until  new  data  supports 
Figure  2. 1.2. 3,  this  result  will  not  be  used  in  this 


research. 

2 . 2  Aerosol  Behavior  Equation 

The  dynamics  of  dispersed,  aerosol  systems  such  as  the 
one  following  a  HCDA  is  described  by  the  aerosol  behavior 
equation.  Derivation  of  this  equation  is  straightforward 
and  is  given  by  Hidy  and  Brock^,  who  also  discuss  its 
solution. 

The  number  of  particles  of  a  given  size  in  the  contain¬ 
ment  atmosphere  at  any  time  can  be  defined  in  terms  of  the 
aerosol  density  distribution  function: 

n(x,t)dx  -  Number  of  aerosol  particles  of  size  x  in 
dx  at  time  t  per  unit  volume  (#/cm3) . 
Applying  this  definition  the  governing  aerosol  rate  equation 
is  written  in  the  following  form: 
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!=■  n(x,t)  =  i  /  dx ' $ (x-x '  ,x)n(x-x '  ,t)n(x',t) 


-  n(x ,  t)  /  dx'4>(x'  ,x)n(x'  ,t) 


-  n(x,t)R(x)  +  S (x , t) 


(2.2.1) 


where 


volume  of  aerosol  particle  with  volume 
equivalent  radius  r^  v 


|  ”  rl,v  (cn,5) 


volume  of  aerosol  particle  with  volume 


equivalent  radius  r2  v 


'  f  "  r2,v  (c”5) 


$(x',x)  ~  Normalized  collision  kernel  predicting 

the  probability  of  collision  between  two 
particles  of  sizes  x'  and  x  (cm3) 


R(x)  ~  Removal  rate  of  particles  of  size  x  due  to 
settling  to  the  containment  floor,  wall 
plating,  leakage,  and  so  forth  (#/sec) 


S(x,t)dx  -  Source  rate  of  particles  of  size  x  in  dx  at 
time  t  per  unit  volume  per  unit  time 
( #/cm^- sec) . 

The  concept  of  volume  equivalent  radius  is  explained  in 
detail  in  Chapter  IV.  For  our  purpose  here,  it  is  the  radius 
of  a  sphere  having  the  same  volume  as  the  irregularly  shaped 
particle . 

The  first  integral  in  Equation  2.2.1  represents  the  for¬ 
mation  rate  of  particles  of  size  x  due  to  collisions  between 
particles  of  size  x'  and  x-x'.  The  integral  is  multiplied  by 
1/2  so  as  not  to  count  each  collision  twice.  The  second 
integral  represents  the  disappearance  rate  of  particles  in  the 
size  range  between  x  and  x  +  dx  due  to  collisions  with  all 
other  particles. 

The  third  term  of  Equation  1.1  is  a  removal  factor  which 
accounts  for  a  decrease  in  n(x,t)  resulting  from  non-colli- 
sional  losses  of  particles  of  size  x.  Losses  include  leakage 
from  the  containment  building,  deposition  on  internal  sur¬ 
faces  due  to  gravitational  settling,  thermophoresis,  and 
Brownian  diffusion,  and  mechanical  removal  by  a  recirculation 
cleanup  system.  The  fourth  term  accounts  for  all  source 
particles  of  size  x. 

Understanding  the  nature  of  the  various  quantities 
described  in  Equation  2.2.1  is  of  considerable  interest.  If 
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should  be  recognized,  however,  that  this  equation  is  not 

strictly  rigorous,  as  it  omits  a  number  of  variables  which 

14 

affect  actual  aerosol  behavior 

To  what  extent  these  variables  will  influence  results 
depend  on  validation  experiments  and  development  of  a  bench¬ 
mark  code.  Work  has  begun  on  HAARM- 3  code  verification 
procedure1 ^ . 

A  benchmark  code  must  describe  all  phenomonological 
mechanisms  believed  to  be  important  in  predicting  the  aero¬ 
sol  density  distribution.  The  resultant  code  is  then  to  act 
as  a  reference  against  which  HAARM- 3  and  other  approximate 
but  more  practical  codes  could  be  tested.  Such  a  benchmark 
code,  called  CRAB,  is  under  development  at  Battelle  Columbus 
Laboratories16 . 

2 . 3  Aerosol  Particle  Coagulation 

Coagulation  is  used  to  describe  the  process  of  growth 
of  aerosol  particles  through  collision  and  their  coalescence 
with  one  another.  Growth  of  aerosol  particles  is  described 
by  such  interaction  mechanisms  as  Brownian  coagulation, 
gravitational  coagulation,  and  turbulent  coagulation.  Since 
relatively  little  is  known  about  the  efficiency  of  coales¬ 
cence,  the  coagulation  efficiency  is  assumed  equal  to  the 
collision  efficiency,  i.e.,  the  probability  of  sticking  after 
collision  is  unity. 
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The  collision  kernel  in  Equation  2.2.1  gives  the  prob¬ 
ability  of  collision  between  particles  of  size  x,x'.  The 
three  dominant  mechanisms  that  make  up  the  collision  kernel 
are  assumed  to  be  separable  and  additive.  The  kernel  can 
be  written: 

<Kx',x)  =  B(x',x)  +  G  (  x  *  ,  x )  +  T(x  '  ,  x)  (2.3.1) 


where 

B(x',x)  ~  Brownian  coagulation 
G(x' ,x)  ~  Gravitational  coagulation 
T(x’  ,x)  -  Turbulent  coagulation. 

The  validity  of  the  assumption  that  the  three  coagula¬ 
tion  mechanisms  are  separative  and  additive  has  not  been 
established.  Synergistic  effects  probably  exist  for  certain 
cases,  but  before  any  of  these  cases  can  be  investigated, 
each  of  the  three  coagulation  processes  must  be  accurately 
defined.  The  functional  form  of  $(x',x)  depends  on  the 
dominant  mechanism  of  coagulation  and  the  Knudsen  number. 
Detailed  discussions  of  the  expressions  for  <j>(x',x)  for  the 
different  coagulation  mechanisms  are  given  in  Hidy  and 
Brock1^.  For  conciseness  only  the  gravitational  coagulation 
term,  G(x',x),  will  be  developed  here  since  this  is  the 
dominant  term  in  the  problem  being  considered  in  this  dis¬ 


sertation. 
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2 . 4  Gravitational  Coagulation 

The  gravitational  coagulation  term,  G(x',x),  of  liquation 
2.3.1  can  be  written  as: 

G(x',x)  =  °G(rl,v,r2,v)  iV(x')  '  V(-X)l  (2.4.1) 


where : 


r.  ~  volume  equivalent  radius  of  aerosol 

*  f  * 


particle  of  volume  x 

r~  -  volume  equivalent  radius  of  aerosol 

L  ,  V 

particle  of  volume  x' 

V(x)  ~  terminal  settling  velocity  of  aerosol 
particle  x 

V(x')  ~  terminal  settling  velocity  of  aerosol 
particle  x' 

°G*-rl  v’r2  v^  ~  gravitational  collision  cross  section 
between  particles  of  radii  r^  and  r2 
undergoing  gravitational  motion. 

The  gravitational  collision  cross  section  defined  in 
Equation  2.4.1  can  be  written  as: 


°G^rl,v,r2,v)  =  7Tyc^rl,v’r2,v^) 


'  "R12  c(rl,v-r2,v> 


(2.4.2) 


(2.4.3) 
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where  y  c  C  r  v»r>  v)  is  the  maximum  horizontal  separation 
between  the  particle  of  size  x  and  x'  that  would  lead  to  a 
collision.  The  term,  R^.  t^ie  sum  particle  radii 

and  e,  the  gravitational  collision  efficiency,  is  defined  by: 


G(rl  V’r^  v>  =  - )2 

A,v  l,v  2,v 


(2.4.4) 


The  gravitational  collision  efficiency,  defined  in 
Equation  2.4.4  is  a  nondimensional  quantity  defining  the 
collision  process  as  a  function  of  the  particle  equivalent 
radii  and  initial  horizontal  separation.  A  gravitational 
efficiency  of  one  means  that  the  center  of  the  smaller 
particle  can  lie  anywhere  in  a  circle  of  radius  r.  +  r? 
and  there  will  be  a  collision  between  the  two  particles.  If 
there  is  hydrodynamic  interactions  between  the  two  particles, 
the  particles  will  have  horizontal  velocity  components  which 
will  tend  to  push  them  apart  and  cause  the  gravitational  col¬ 
lision  efficiency  to  be  less  than  one.  If  the  efficiency  is 
equal  to  zero,  the  smaller  aerosol  particle  would  have  to  be 
initially  positioned  directly  under  the  larger  aerosol  parti¬ 
cle  if  any  collision  is  to  occur,  and  even  then  there  is  the 
possibility  that  the  particles  will  miss  each  other.  In  the 
actual  case,  the  gravitational  collision  efficiency  will 
depend  on  the  hydrodynamic  interactions  which  will  be  a 
function  of  such  quantities  as  aerosol  particle  size,  density, 
and  shape  and  contaminant  atmospheric  properties. 


2.5  Problem  Statement 


The  purpose  of  this  research  is  to  investigate  the 
gravitational  collision  efficiency  for  post-HCDA  LMFBR  aero¬ 
sol  particles.  As  discussed  in  Section  2.1,  agglomerates  of 
sodium  oxides,  plutonium  and  uranium  oxides  are  branched 
chain-like  structures.  To  date  collision  efficiency  calcu¬ 
lations  are  based  on  the  assumption  that  particles  are 
spherically  shaped.  Results  from  a  sensitivity  analysis 
of  the  HAARM-3  Code  have  shown  that  this  model  is  sensitive 
to  the  gravitational  collision  efficiency,  density  correc¬ 
tion  term,  and  two  particle  shape  factors  designed  to  correct 
for  the  deviation  from  sphericity  of  the  agglomerates. 

The  theoretical  determination  of  the  gravitational  col¬ 
lision  efficiency  requires  solving  the  full,  steady  state 
Navier- Stokes  equation  of  motion.  Given  the  velocity  fields 
and  drag  coefficients  of  the  particles  involved,  the  super¬ 
position  method  can  be  used  to  approximate  the  actual  hydro- 
dynamic  interaction  between  the  two  particles.  Numerical 
integration  methods  are  then  used  to  solve  the  first-order, 
coupled  differential  equations  which  describe  the  particles' 
trajectories.  The  critical  trajectory  which  results  in  a 
grazing  collision  gives  y  .  At  present,  the  theoretical  work 
on  collisional  mechanism  and  efficiency  appears  more  advanced 
than  the  experimental  information.  This  is  the  state  of 
affairs  for  the  work  in  this  dissertation  also. 
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The  contributions  of  this  research  may  be  summarized 
as  follows: 

1.  Investigation  of  the  nonspherical  gravitational 
collision  kernel.  To  date  all  aerosol  behavior  codes  have 
used  very  simple  collision  kernels  based  on  the  assumption 
that  LMFBR  aerosols  were  spherical  in  shape.  Any  correction 
for  nonspherical  particles  was  based  on  single  particle 
shape  factors. 

2.  Investigation  of  the  slow  convergence  rate  of  the 
Navier-Stokes  equation  as  a  function  of  Reynolds  number  and 
axis  ratio. 

3.  Development  of  a  generalized  computer  code,  NGCEFF, 
written  to  solve  the  Navier-Stokes  equation,  the  equations 
of  aerosol  motion  and  to  determine  the  matrix  parameters 
affecting  convergence. 

4.  Development  of  several  function  programs  using  the 
collision  efficiencies  determined  from  the  spherical  colli¬ 
sion  efficiency  program,  GCEFF. 

5.  Formulation  of  function  routines  for  nonspherical 
particles  which  will  facilitate  the  use  of  the  results  of 
NGCEFF. 

6.  A  preliminary  investigation  of  the  effects  of  slip 
on  collisions  efficiency  when  kinetic  corrections  are  made 
to  the  Stokes  approximation. 


CHAPTER  III 


REVIEW  OF  PREVIOUS  WORK 


3.0  Introduction 

The  preceding  chapters  introduced  the  research  problem 
associated  with  the  calculation  of  the  gravitational  colli¬ 
sion  efficiency.  Most  of  the  interest  i:;  this  quantity  in 
the  past  has  been  in  the  atmospheric  sciences  in  conjunction 
with  the  prediction  of  the  growth  rate  of  raindrops  in  the 

clouds.  One  notable  exception  is  the  recent  work  of  Pertmer 
1 8 

and  Loyalka  ,  who  investigated  the  gravitational  collision 
efficiency  of  aerosols  formed  following  an  HCDA.  In  comput¬ 
ing  the  collision  kernels,  researchers  in  the  atmospheric  and 
nuclear  sciences  relied  on  the  work  of  Stokes,  Oseen,  Stimson 
and  Jeffery,  Hocking  and  others  to  determine  the  flow  field 
around  a  single  sphere  and  the  drag  forces  between  the  hydro- 
dynamic  interacting  spheres.  Consequently,  much  is  known 
about  the  collision  efficiencies  of  two  hydrodynamically 
interacting  rigid  spheres,  while  very  little  is  known  about 
the  collision  efficiencies  of  nonspherical  aerosols. 

In  this  chapter,  the  previous  investigations  of  the  flow 
field  and  drag  calculations  for  oblate  spheroids,  and  the 
collision  efficiency  for  nonspherical  aerosols  are  reviewed. 
Appropriate  references  will  be  made  to  the  contributions  of 
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those  researchers  who  worked  with  spherical  particles  but 
for  a  detailed  review  of  this  body  of  literature,  the  reader 
should  consult  other  references^  ^  . 

3. 1  Flow  Field  Approximation 

To  determine  the  gravitational  collision  kernel,  the 
velocity  fields  and  drag  forces  between  two  interacting  aero¬ 
sol  particles  must  be  known.  As  stated  in  previous  chapters, 
an  oblate  spheroid  will  be  used  to  model  the  nonspherical 
particle. 

In  principle,  the  velocities  and  forces  on  nonspherical 
aerosols  can  be  determined  by  solving  the  continuity  and 
Navier-Stokes  equation  of  momentum.  Since  these  partial  dif¬ 
ferential  equations  are  nonlinear,  closed  form  solutions  are 
available  for  only  special  cases. 

The  only  exact  solutions  of  a  two  particle  problem  is 

4  9 

that  of  Stimson  and  Jeffery  ,  for  the  Stokes  motion  of  two 
spheres  parallel  to  their  line  of  centers.  The  system  of 
bipolar  coordinates  employed  is  unique  in  that  it  permits 
one  to  satisfy  simultaneously  the  boundary  conditions  on  two 
external  spheres. 

For  the  case  of  a  sphere  colliding  with  a  nonsperhical 
body,  it  is  not  generally  possible  to  find  coordinate  systems 
which  permit  simultaneous  satisfaction  in  all  the  boundary 
conditions.  Accordingly,  numerical  schemes  must  be  sought 
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whereby  the  boundary  value  problem  may  be  solved  by  consider¬ 
ing  boundary  conditions  associated  wi'th  one  particle  at  a 
time.  Even  numerical  schemes  are  limited  in  their  ability 
to  predict  flow  conditions  around  nonspherical  particles. 
Consequently,  nonspherical  particles  must  be  chosen  which 
permit  investigation  of  viscous  flow  around  objects  of  ro¬ 
tational  symmetry.  An  oblate  spheroid  possesses  such  rota¬ 
tional  symmetry. 

3.1.1  Steady  Flow  Past  Oblate  Spheroids 

The  problem  of  calculating  numerically  the  velocity 

field  and  drag  coefficient  for  viscous  fluid  flowing  past  a 

single  sphere  is  well  documented  in  the  literature. 

20 

Jensen  ,  using  the  Stokes  stream  function  and  the  concept 

of  vorticity,  successfully  applied  the  finite  difference 

method  to  solve  the  Navier- Stokes  equation  for  flow  past  a 

sphere  at  Reynolds  numbers  up  to  40.  Later,  Rimon  and 
21 

Cheng  obtained  the  transient  uniform  flow  around  a  sphere 

for  Reynolds  numbers  to  1000.  However,  their  results  for 

the  surface  pressure  distribution  at  steady  state  conditions 

are  a  variance  with  those  of  Jensen  and  of  Pearcey  and 
2  2  2Z 

McHugh  .  Rimon  and  Lugt  ,  using  similar  methods  as  those 
of  Rimon  and  Cheng,  arrived  at  the  steady  state  solution  for 
flow  past  oblate  spheroids  at  Reynolds  numbers  of  10  and  100. 
They  showed  that  the  earlier  work  of  Michael^  does  not 
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correctly  treat  the  sharp  edge  of  the  disk  for  axis  ratios 
0.2  and  0.5,  where  the  axis  ratio  is  defined  as  the  minor 
axis  divided  by  the  major  axis.  However,  the  works  of 

2  S 

Rimon  and  Lugt  are  also  suspect,  as  discussed  by  Pruppacher 
et  al . ,  due  to  inherent  difficulties  in  reaching  steady-state 
conditions  by  solving  the  time-dependent  Navier- Stokes  equa¬ 
tion  of  motion  and  problems  associated  with  their  numerical 
method  of  solution,  which  also  plagued  their  treatment  of 
the  spherical  case. 

For  oblate  spheroids,  the  best  results  are  those  reported 

7  2  7 

by  Pitter  et  al .  ,  and  Masliyah  and  Epstein  .  Pitter 

performed  calculations  for  spheroids  of  axis  ratio  0.05  and 
0.2  and  Reynolds  numbers  between  0.1  and  100  while  Masliyah 
and  Epstein's  results  were  for  axis  ratio  0.2  and  0.5  and 
Reynolds  numbers  between  1  and  100.  Figures  3. 1.1.1  and 
3. 1.1. 2  show  the  influence  of  axis  ratio  and  Reynolds  num¬ 
bers  on  the  flow  around  oblate  spheroids.  The  vorticity  is 
generated  upstream  and  is  carried  by  the  fluid  around  the 
spheroid  to  considerable  distances  downstream.  It  is  inter¬ 
esting  to  note  that  the  streamlines  of  magnitude  4.0  tend  to 
be  less  curved  with  increasing  Reynolds  number,  indicating 
that  the  flow  becomes  undisturbed  at  a  shortened  distance 
from  the  spheroid  as  Reynolds  number  increases.  This  result 
correlates  with  the  results  that  Pitter  obtained  when  he 
calculated  the  collision  efficiency  between  an  oblate 
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Figure  3. 1.1.1  Streamlines  for  Oblate  Spheroids 
with  Axis  Ratio  0.2  Vorticity  Lines  Are 
Indicated  by  Dash  Lines  Reynolds  Number 
(a)  10,  (b)  50,  (c)  10027 
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spheroid  with  an  axis  ratio  0.05  and  a  colliding  sphere  and 
found  that  collision  occurred  on  the  outer  edges  of  the 
spheroid  but  not  on  the  interior  of  the  spheroid  for  higher 
Reynolds  numbers.  This  surprising  result  will  be  discussed 
later  in  this  chapter. 

Equally  important  were  the  findings  of  Pitter  for  the 

drag  on  oblate  spheroids.  He  determined  that  as  the  Reynolds 

number  becomes  less  than  one,  the  drag  is  approximated  better 

by  the  Oseen  drag  expression  for  oblate  spheroids  than  by  way 

of  creeping  flow  drag  expression.  This  is  consistent  with 

2  8 

the  finding  of  LeClair  et  al .  for  spheres,  and  with  the 

2  5 

findings  of  Pruppacher  et  al .  for  cylinders.  Furthermore 
for  Reynolds  numbers  less  than  0.01  Pitter  claims  that  the 
drag  on  a  thin  oblate  spheroid  is  given  accurately  by  the 
Oseen  drag  expression. 

Related  to  these  findings  of  Pitter,  is  an  important 

29 

result  presented  by  Klett  for  two  particle  motion.  He 
found  that  strong  relative  motions  exist  between  two  parti¬ 
cles  at  very  small  Reynolds  numbers,  where  according  to 
Stokes  theory  the  viscous  forces  should  prevent  this  motion. 
This  indicates  that  the  inertia  terms  of  the  full  Navier- 

Stokes  equation  cannot  be  neglected  and  possibly  why  Pertmer 
1  8 

and  Loyalka  found  that  the  Oseen  approximation  agreed 

4  8 

better  with  the  experimental  data  of  Tu  and  Shaw  for 
cloud  droplets  than  the  results  of  other  researchers  who 
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used  Stokes’  approximation.  Another  possibility  is  dis¬ 
cussed  in  Appendix  A3,  Kinetic  Corrections  to  Aerosol  Gravi¬ 
tational  Collision  Efficiency. 

Nevertheless,  there  is  enough  doubt  about  Stokes'  ap¬ 
proximation  for  calculations  related  to  the  gravitational 
collision  efficiency  that  any  solution  should  use  the  full 
Navier-Stokes  equation.  This  will  be  discussed  further  when 
the  so-called  superposition  method  is  reviewed. 


3.1.2  Steady  Flow  Past  Spheres 


This  topic  has  been  investigated  both  analytically  and 

numerically.  Analytical  solutions  are  limited  to  Stokes, 

Oseen  and  variations  of  them.  For  higher  Reynolds  numbers, 

empirical  and  numerical  methods  must  be  used.  For  their 

2  7 

research,  Masliyah  and  Epstein  showed  that  an  oblate 
spheroid  with  an  axis  ratio  0.999  can  be  considered  to  be  a 
sphere  for  Reynolds  numbers  1  to  100.  This  is  one  advantage 
of  using  oblate  spheroids  to  model  nonspherical  particles. 


3.2  The  Superposition  Model 


In  order  to  model  the  hydrodynamic  interaction  of  two 

29 

spheres  with  as  much  physical  realism  as  possible,  Klett  , 

1 8 

Pertmer  and  Loyalka  ,  and  others  used  analytical  methods, 
based  on  Stokes  and  Oseen  equations,  which  simultaneously 
satisfied  the  changing  boundary  conditions  for  interacting 
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spheres.  While  such  analytical  treatment  is  highly  desirable 
for  the  problem  being  considered  here,  the  complex  geometry 
associated  with  the  interaction  between  an  oblate  spheroid 

and  sphere  presently  prevents  it. 

30  31  32 

Langmuir  ,  Shafrir  and  Neiburger  ,  Pearcey  and  Hill  , 

3  3  18 

Shafrir  and  Gal-Chen  ,  and  Pertmer  and  Loyalka  have  used 
the  so-called  superposition  model  to  compute  the  collision 
efficiencies  of  aerosols.  In  this  model,  the  motion  of  each 
particle  is  superpositioned  onto  a  fluid  velocity  field  and 
the  drag  force  is  assumed  to  be  taken  into  account  by  the 
velocity  difference  between  the  particle's  velocity  and  the 
fluid  velocity  field  of  the  second  particle  superpositioned 
onto  the  first  particle. 

In  order  to  utilize  this  model,  the  flow  fields  of  each 

1 8 

particle  must  be  defined.  Pertmer  and  Loyalka  employed 
Stokes'  approximation  to  the  Navier- Stokes  equation  to  obtain 
the  flow  field  around  one  aerosol  particle  as  it  falls  in 
isolation.  This  was  one  of  the  options  available  in  their 
code,  GCEFF.  They  reported  that  the  superposition  method 
and  the  analytical  model  using  Stokes'  approximation  gave 
similar  collision  efficiencies.  This  seems  to  indicate  that 
the  superposition  method  is  a  good  approximation  to  be  used 
to  estimate  the  drag  force  terms,  if  certain  conditions  are 
satisfied.  These  are  discussed  below. 


3  2 

Pearcey  and  Hill  utilized  the  Oseen  approximation  in 

31 

order  to  determine  the  flow  field.  Shafrir  and  Neiburger  , 

on  the  other  hand,  solved  the  Navier- Stokes  equation  to  eval- 
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uate  fluid  flow  velocities,  as  did  Lin  and  Lee  .  A  compar¬ 
ison  of  all  of  these  results  is  possible  but  since  the 
resulting  equations  of  motion,  necessary  to  calculate  the 
gravitational  collision  efficiency,  possess  the  problem  of 
mathematical  stiffness,  it  is  difficult  to  isolate  differ¬ 
ences  due  to  drag  force  terms  and  differential  equation  inte¬ 
grator  routines  used  by  the  various  investigators. 

7  C 

Pitter  and  Pruppacher  used  the  superposition  method 
to  calculate  collision  efficiencies  between  falling  ice 
crystals  and  raindrops.  They  developed  an  extensive  argument 
which  justified  the  use  of  this  method.  Basically  they  found 
that  a  reasonable  quantity  for  evaluating  the  effect  of  the 
smaller  particle  flow  field  on  the  larger  particle  is  the 
mass  ratio  of  the  two  particles.  For  interactions  investi¬ 
gated  by  them  this  quantity  never  exceeded  0.05  even  though 
the  velocity  ratio  of  interacting  water  drop  and  ice  crystal 
was  as  high  as  0.8.  This  is  important  since,  when  bodies 
have  similar  fall  velocities,  the  results  are  questionable 
due  to  the  length  of  time  for  interactions  to  take  place. 
Another  important  reason  to  keep  the  mass  ratio  small  is  to 
prevent  tilting  of  the  crystal. 
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Thus  the  possibilities  available  to  represent  the  flow 
field  include  Stokes  and  Oseen  approximation  and  solutions 
to  the  full  Navier- Stokes  equation.  As  discussed  earlier, 
Stokes'  solution  for  flow  around  oblate  spheroids  is  ques¬ 
tionable.  Oseen's  solution  is  also  limited  to  low  Reynolds 
number.  It  was  decided  that  numerical  solutions  to  the  full 
Navier-Stokes  equation  offered  the  greatest  capability. 


Gravitational  Collision  Efficiencies  for  Nonspherical 


Particles 


If  the  problem  of  particle  deposition  of  fibrous  filters 
is  omitted,  then  it  has  been  only  recently  that  some  attempts 
to  obtain  collision  efficiencies  for  nonspherical  bodies  has 
appeared  in  the  literature.  The  most  recent  is  by  Pitter  and 
Pruppacher^5  who  numerically  investigated  the  hydrodynamic 
interaction  between  simple  ice  plate,  modeled  as  oblate 
spheroids  of  axis  ratio  0.5  and  spherical  water  drops.  The 
spheroids  had  semi-major  axis  lengths  between  147  and  404  ym 
and  drops  with  radii  up  to  53  ym.  As  indicated  earlier, 
their  investigation  revealed  the  surprising  result  that  water 
spheres  of  certain  sizes  were  found  which  would  not  collide 
inside  some  inner  critical  trajectory  offset,  termed  y.  .  , 
which  is  less  than  the  initial  offset  of  the  critical  grazing 
trajectory,  denoted  as  y  .  In  other  words,  if  ym^n  *  0, 
the  collision  domain  beneath  the  spheroid  can  be  pictured 


as  a  circular  cylinder  with  its  axis  of  symmetry  of  the 
spheroid.  However,  when  ymin  >  0,  the  collision  domain  can 
be  pictured  as  a  circular  annulus  of  radii  y  .  and  y  . 

Figure  3.3.1  shows  the  variation  of  yc  and  ymjn  for  an 
ice  plate  with  a  major  axis  equal  to  160  pm  for  different 
sphere  radii.  Figure  3.3.2  displays  collision  efficiencies 
of  several  ice  plates  interacting  with  water  droplets.  The 
interesting  feature  of  Figure  3.3.2  is  the  cutoffs  at  large 
and  small  raindrops  sizes.  The  decrease  in  droplet  size, 
which  can  collide  with  increasing  ice  plate  size  is  due  to 
the  inertia  of  the  small  droplet  moving  it  across  streamlines 
and  larger  terminal  velocities  for  larger  spheroids.  This 
effect  can  also  be  deduced  from  Figure  3. 1.1.1,  as  can  the 
possibility  of  having  an  annular  collision  domain.  Basically 
a  sphere  tends  to  fall  at  its  terminal  velocity  with  respect 
to  the  fluid  immediately  surrounding  it,  at  least  this  is  the 
case  for  most  collision  efficiency  models.  Also,  a  sphere's 
acceleration  is  proportional  to  the  difference  between  the 
velocity  the  sphere  tends  to  reach  and  its  actual  velocity. 

If  for  y  >  yc  and  the  water  sphere  is  able  to  accelerate  to 
the  terminal  velocity  of  the  oblate  spheroid  of  ice  before 
they  collide,  the  horizontal  forces  will  move  the  drop  around 
the  ice  plate  before  a  collision  occurs.  Naturally  the 
largest  y  at  which  this  occurs  is  called  ym^n» 


Figure  3.3.2  Collision  Efficiency  as  a  Function  of 
Oblate  Spheroid  Size,  a^,  and  Sphere  Radius,  as. 
Oblate  Spheroid  Sizes:  (1)  160,  (2)  194, 

(3)  213,  (4)  289,  (5)  404  micrometers. 


The  cutoffs  at  larger  water  droplets  sizes  is  more 
difficult  to  explain.  Pitter  and  Pruppacher  gave  no  explan¬ 
ation  for  this  effect  but  it  may  be  due  to  the  increased 
interaction  time  which  would  allow  the  horizontal  forces 
to  act  longer  on  the  raindrop.  Another  possibility  may  be 
the  integration  routine  used  to  solve  the  system  of  differ¬ 
ential  equations,  which  are  stiff  as  showed  by  Pertmer  and 
Loyalka*® . 


CHAPTER  IV 


MATHEMATICS  OF  NONSPERICAL  AEROSOL  GRAVITATIONAL  COAGULATION 

4 . 0  Introduction 

The  gravitational  collision  efficiency  has  been  dis¬ 
cussed  and  its  importance  in  the  calculation  of  the  aerosol 

density  distribution  function  has  been  established  by  sen- 

1 7 

sitivity  analysis  .  Many  researchers  feel  that  the  numeri¬ 
cal  methods  necessary  to  solve  the  aerosol  rate  equation 
are  sufficiently  well  advanced  and  what  really  is  needed  at 
present  are  fairly  accurate  collision  kernels  for  the  aero¬ 
sol  rate  equation. 

In  this  chapter,  the  mathematics  describing  the  colli¬ 
sion  process  between  an  oblate  spheroid  and  a  sphere  is  pre¬ 
sented.  Where  possible,  generalized  equations  are  derived 
which  can  apply  to  any  nonspherical  particle  or  particles. 
However,  to  solve  the  stated  research  problem,  generalized 
equations  must  be  abandoned  for  the  detailed,  specific 
equations.  It  is  hoped  that  the  concepts  and  technique 
presented  can  be  used  for  other  shaped  particles. 

4 . 1  Gravitational  Coagulation  and  Shape  Factors 

As  defined  in  Chapter  2,  Equation  2.4.1,  the  gravita¬ 
tional  collision  kernel  of  two  aerosol  particles  of  volume 
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(4.1.1) 


(4.1.2) 

where  e^fr^  v»r2  v)  t^e  nonsP^erical  gravitational  colli¬ 
sion  efficiency  for  particles  x^  and  x2  with  terminal  set¬ 
tling  velocities  VCx-^)  and  V(x2)  . 

Mass  equivalent  radii  are  defined  by  the  relationship: 

ri,m  =  (3mi/4TTP.)1/3  i  -  1,  2  (4.1.3) 

where  m^  and  p^  are  the  mass  and  material  density  of  the 

i**1  particle.  In  the  case  of  agglomerates,  the  actual 
(measured)  agglomerate  density  pT,  is  related  to  the  material 
density  by: 

pi  =  aipi  i  =  1,  2  (4.1.4) 

f  h 

where  is  called  the  density  correction  factor  for  the  iv 
particle.  The  volume  equivalent  radii  are  then  related  to 
the  mass  equivalent  radii  by  the  relationship: 


equivalent  radii,  r.  and  r?  ,  is  written  as: 

1  >  V  L  j  V 

G(x^,x2)  =  °q  ( r2.  t  v »  r2  ,  v^  l^(xi^  '  ^  (x2)  I 


or  in  its  expanded  form: 


G(xlfx2)  =  TT(r-,  v  +  r,  J  v,r?  J 


1 ,  v  ‘2,vJ  N'-  1  ,v*  2  ,vJ 


I  V ( x , )  -  V(x7) 
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ri,v 


a/1/3  r • 


(4.1.5) 


i  =  1,  2 


r . 
i,v 


f  r . 


(4.1.6) 


where  Equation  4.1.6  defines  a  "shape  factor"  referred  to  as 
the  effective  diameter  factor.  For  the  remainder  of  this 
chapter,  the  notation  of  the  it^1  particle  will  refer  to  par¬ 
ticle  Xj  or  x 2 . 

There  are  several  other  "shape  factors"  which  are  quoted 
in  the  literature  and  applied  in  aerosol  codes  such  as 
HAARM-32.  Since  particles  encountered  in  the  environment  are 
generally  irregular  in  shape,  the  diameter  assigned  to  an 
aerosol  particle  usually  depends  on  the  manner  in  which  it  is 
measured.  A  fundamental  concept,  called  sphericity,  measures 
the  deviation  of  an  individual  particle's  shape  from  that 
characterized  by  a  sphere.  It  is  defined  as  the  ratio  of  the 
surface  area  of  a  sphere  of  the  same  volume  as  the  particle 
to  the  actual  surface  area,  F,  of  the  particle.  Mathemati¬ 
cally,  sphericity,  $,  is  expressed  as: 


4>  = 


4 7i  (rv)  2 


(4.1.7) 


where  ry  is  the  volume  equivalent  radius  of  a  sphere  of  the 
same  volume  as  the  original  particle. 


42 


For  investigations  of  the  motion  of  aerosols  the  well 
known  Stokes  law  is  often  used.  Stokes  showed  that  the 
resistance  to  the  motion  of  a  spherical  particle  through  a 
viscous  fluid  will  be  proportional  to  the  product  of  the 
viscosity  of  the  fluid  and  the  velocity  and  diameter  of  the 
particle,  provided  the  fluid  medium  is  infinite  in  extent  and 
that  the  inertia  forces  are  negligible  compared  to  the  vis¬ 
cous  forces.  If  the  shape  of  a  particle  is  not  spherical, 
the  single  dimension  of  diameter  will  no  longer  be  sufficient 
to  describe  the  resistance  of  the  fluid  to  the  motion  of  this 
particle.  An  aerodynamic  diameter  d&  has  been  utilized  in 
this  case  and  is  defined  as  the  diameter  of  a  spherical  par¬ 
ticle  with  unit  density,  pQ,  with  the  same  settling  velocity 
as  that  measured  for  a  given  nonspherical  particle.  This 
diameter  tends  to  obscure  shape  effects  rather  than  explain 
them,  since  often  neither  the  shape  nor  the  size  of  a  parti¬ 
cle  is  measured.  Mathematically,  aerodynamic  radius,  r. 

1  9  a 

is  given  by  the  formula: 


Vi.a 


2  po«  ri , a 
9  V  D(ri>a) 


(4.1.8) 


where  V.  is  the  actual  (measured)  settling  velocity  for 
i  ,  a 

f.  L 

the  i  particle,  g  is  the  gravitational  constant,  \i  is  the 

fluid  viscosity,  and  D(r.  _)  is  the  corrected  drag  coeffi- 

i  ,  a 

cient  defined  by: 
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Re 


DCri,a)  =  ^ri,a)  ff 


(4.1.9) 


where  <Kr.  _)  is  the  measured  drag  coefficient  and  Re  is  the 
i ,  a 

Reynolds  number.  For  aerosols,  the  variable  D(r.  )  is 

l ,  a 

approximated  by: 


D(ri,a)  = 


1.0  +  0.13  Re.  0-85 

i  .a 


C(ri,a> 


(4.1.10) 


where  C(r.  )  is  the  Cunningham- Stokes  correction  term, 
i ,  a 

applicable  when  the  Knudsen  number  is  not  negligible. 

Since  the  aerodynamic  radius  tends  to  obscure  shape 
effects,  another  shape  factor  based  on  a  ratio  of  settling 
velocities  has  been  defined.  If  the  volume  equivalent  radius 
and  actual  agglomerate  density,  defined  by  Equations  4.1.5 
and  4.1.4,  respectively,  are  used  in  Stokes  law,  a  settling 
velocity  different  from  the  measured  settling  velocity  would 
be  calculated.  In  this  context,  the  volume  equivalent  rad¬ 
ius  is  sometimes  called  the  Stokes  radius  in  order  to  con¬ 
form  with  its  use.  The  Stokes  settling  velocity,  V.  ,  based 
on  the  volume  equivalent  radius  is  given  by: 


Vi,s 


2  pi§ 


r? 


i  ,v 


D(rijV) 


(4.1.11) 
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where  D(r.  )  is  the  correction  drag  coefficient  for  the 

1  f  * 

Stokes  radius.  At  this  point  the  aerodynamic  or  dynamic 
shape  factor  ,  k^,  is  defined  as  the  ratio  of  the  Stokes  set¬ 
tling  velocity  to  the  aerodynamic  (actual  or  measured)  set¬ 
tling  velocity.  Thus: 


< . 
1 


Vi,s  =  /i.V) 2  U(~ri.a:) 

vi,a  po  ri,a  D(ri>v) 


(4.1.12) 


or  applying  Equation  4.1.10 


K 


i 


(lizV) 

ri,a 


C(ri>v)  (1  *  0.13  Re^;85) 
C(ri,a}  t1  +  °-13  Rei%85) 


(4.1.13) 


Equation  4.1.13  points  out  that  the  drag  correction  factor 

t*  H 

depends  on  the  Reynolds  number  for  the  i  particle  and  on 
the  radius  used  to  characterize  the  particle.  In  the  major¬ 
ity  of  cases,  the  approximate  shape  factor,  k^,  defined  by: 


K  . 


1 


p  ■  r  - 
=  li  (-±tX) 

po  ri,a 


2 


C^i,a) 
1  cCri,v) 


(1  +  0.13  Re.0-85) 
_ _  ~  >  _ _ 

(1  +  0.13  Re-0;85 

i ,  a 


(4.1.14) 


can  be  used. 
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# 


G(x1,x2)  -  ^(r1>m  +  r2>m)  eN(ri>m»r2>mJ 


1/3 

2  g  tti  Pi 
9  u  <i 


2  cvAai/3 

rl.«  “  *  ^rT1-) 

1  ,m 


1/3  P  ^  K-j  ^ 

<4>  (pf>  ^ 


V“2/3 

(1  ♦  -^--Z--) 

r2,m 


(4.1.17) 


The  relationship  between  the  gravitational  collision 
efficiency  based  on  mass  equivalent  radii  and  volume  equiva 
lent  radii  is: 


eN(rl,m’r2,n.> 


a 


TJT 

1 


a.  1/3 

(r,  +  ( — )  •  r0  ) 
^  i  m  '  2  ,m' 


l,m  '•a- 


(rl,m  *  r2,m> 


(4.1.18) 


Since  density  correction  factors  are  less  than  one, 

£N(rl,m’T2,m)  >  eN(rl,v’  ^v5  when  yc  (see  Ecluation  2-4-2) 
is  measured  or  calculated  based  on  the  volume  or  mass  equiv¬ 
alent  radii  of  the  colliding  particles.  In  fact  if  *  a2, 
the  relationship  is  simply: 


'N 


(rl,m,rl,m 


)  = 


2/3 


cN(rl,v'r2,v) • 


(4.1.19) 
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Equations  4.1.16  and  4.1.17  have  done  little  insofar  as 
allowing  the  gravitational  collision  kernel  to  be  evaluated 
since  the  nonspherical  gravitational  collisional  efficiency, 
eN,  is  still  unknown.  If  one  assumes  that  both  particles  of 
size  x,  and  x?  are  spheres  with  radii  r,  „  and  r0  and  that 

A  *•  1  *  V  Z  }  V 

calculations  give  a  maximum  horizontal  separation,  ycs,  which 
determines  a  spherical  gravitational  collision  efficiency, 
es(rl  V’r2  v)  »  t*ie  eS  written  as: 


£S<rl,vr2,v>  ■  <r 


l,v  *2,v 


(4.1.20) 


However,  the  nonspherical  gravitational  collision  effi¬ 
ciency,  eN>  will  in  general  differ  from  Eg  because  the  maxi¬ 
mum  horizontal  separation  will  be  influenced  by  the  geometric 
projection  area  and  hydrodynamic  effects  of  the  nonspherical 
particles.  Equations  4,1.20  is  thus  written  as: 


cN(rl,Vr2,v>  " 


(rl,v  *  r2,v> 


(4.1.21) 


From  Equations  4.1.20  and  4.1.21,  the  gravitational  colli¬ 


sion  shape  factor  8,  is  defined: 

eN(rl,Vr2,v>  '  SEs(rl,Vr2,v> 

8  ■  (I™)2 

ycs 


(4.1.22) 


(4.1.23) 
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Expression  4.1.22  can  now  be  evaluated  by  using  the  spherical 
collision  efficiency  results  of  Pertmer  "  and  a  range  of  B 
values  for  different  LMFBR  aerosol  particles.  Before  dis¬ 
cussing  one  method  to  calculate  8  values,  the  HAARM-3  expres¬ 
sion  for  the  gravitational  collision  kernel  will  be  reviewed. 

Briefly,  HAARM-3  expression  is  derived  from  Equation 
4.1.16  by  implicitly  assuming: 


a  ^  =  a  =  a 


Pi  =  P2  =  p 


<1  =  k2  =  < 


(4.1.24) 


which  reduces  this  expression  to: 


G(x1,x2)  =  71  ( ri  >  v  +  r2,v)  eN  ^  rl ,  v  ’  r2  ,  v^  9  u  kP 


rl  v2(1  +  C  7^—)  -  r2  v2 

i,v  r! ,  v 


(i  +  c 

2,v 


(4.1.25) 


or  by  rearranging  terms: 


G(x1,x2)  =  |  es^rl  ,v,r2,v;> 

'  l(rl,v  +  r2,v)3irl,v  ‘  r2,v' 

+  AC(rl,v  +  r2,v)2  lrl,v  ‘  r2  ,  v I ^  *  (4.1.26) 
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A  similar  expression  can  be  derived  for  the  collision  kernel 
based  on  mass  equivalent  radii.  This  is  easily  done  by  using 
Equation  4.1.17,  remembering  conditions  established  by  Equa¬ 
tion  4.1.24,  and  using  the  definition  of  8,  the  kernel  based 
on  mass  equivalent  radii  becomes: 


G(x^ , X£) 


2_  Tigp  q^'//^8 
9  V  X 


ES(rl,m-r2,m) 


*  r2,m>  lrl,m  '  r2,ml 


♦  a1/3  XC(r.  *  r,  ) 2 
1  ,m  2  ,m' 


(4.1.27) 


At  this  point  it  is  worthwhile  to  compare  the  expression 

? 

used  in  the  HAARM- 3  code  and  the  shape  factors  used  to  cor¬ 
rect  for  nonspherical  particles.  HAARM- 3  defines  a  dynamic 
shape  factor,  x»  as  the  square  of  the  ratio  of  the  volume 
equivalent  diameter  to  the  aerodynamic  diameter  (based  on  p' 
instead  of  p  ) .  Also  a  collision  shape  factor,  y,  is  defined 
as  the  ratio  of  the  so-called  "collision  diameter",  dc,  to 
the  volume  equivalent  diameter,  dy.  It  is  not  clear  whether 
dy  is  based  on  only  one  or  both  particles.  The  HAARM-3 
expression  is: 


i 
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G(x,  ,x7)  =  | 

1’  2  9  UX 


-l/3„2 


ec(r,  ) 

S v  1  ,m’  2  ,nr 


[  (r ,  +  r~  )  |  r.  -  r0 

L  1 ,m  2  ,nr  1  1 ,m  2  ,m 


+  AC(r,  +  r?  )2 

1 ,m  2  ,m' 


(4.1.28) 


The  differences  between  Equations  4.1.27  and  4.1.28  are  due 
to  the  factor  a  and  y  .  The  density  correction  factor  is 
easily  explained  by  the  incorrect  assumption  that 

eSCrl,m»r2,m)  =  £S(rl,v’  r2,v)'  when  U  is  correctly  8iven 
by  Equation  4.1.19.  The  introduction  of  y  is  more  serious 

because  y  lacks  theoretical  basis  in  view  of  the  defining 

relationship  and  the  final  HAARM-3  expression.  It  appears 

that  y  is  used  to  represent  the  same  quantity  as  8  but  is  not 

clear  how  y  would  be  calculated  or  determined  experimentally. 

The  remainder  of  this  chapter  will  provide  insight  into 
the  nature  of  g  for  irregularly  shaped  aerosols  by  first  con¬ 
sidering  the  case  of  oblate  spheroids  of  different  axis 
ratios.  Besides  providing  a  suitable  coordinate  system,  an 
oblate  spheroid  of  varying  axis  ratio  will  provide  a  general 
family  of  possible  aerosol  shapes  by  which  the  collision  pro¬ 
cess  can  be  studied.  It  is  hoped  that  the  methods  developed 
in  this  research  will  provide  the  basis  needed  to  study  more 
complex  shapes. 
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4 . 2  Coordinates  Systems 

As  noted  in  previous  chapters,  coordinate  systems  must 
be  chosen  which  facilitate  the  solution  of  the  Navier-Stokes 
equation.  This  not  only  allows  the  boundary  conditions  to 
be  specified  easily  but  simplifies  the  Navier-Stokes  equation 
expressed  in  the  particular  coordinate  system  chosen.  Oblate 
spheroidal  coordinates  are  the  natural  choice  for  the  transla¬ 
tion  of  any  ellipsoid  parallel  to  a  principal  axis.  Also 
ellipsoids  can  be  varied  from  a  sphere  (axis  ratio  *  1.0)  to  a 
circular  disk  (axis  ratio  =  0.0)  and  thus  a  range  of  nonspher- 
ical  particles  is  available. 

■7 

Happel  and  Brenner  describe  many  orthogonal  curvilinear 
coordinate  systems  associated  with  bodies  of  revolution  in 
terms  of  cylindrical  coordinates.  Since  flow  is  assumed  to  be 
axisymmetrical ,  cylindrical  coordinates  offer  the  advantage 
that  the  interaction  between  nonspherical  bodies  can  be  cal¬ 
culated  using  the  velocity  fields  associated  with  each  parti¬ 
cle  after  being  transformed  from  their  respective  coordinate 
system. 

4.2.1  Oblate  Spheroidal  Coordinates 

Happel  and  Brenner‘S  define  the  oblate  spheroidal  coor¬ 


dinate  system  based  on  the  transformation: 
z  +  ip  =  c  sinh(5  +  in) 


(4. 2. 1.1) 
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where  z  and  p  are  cylindrical  coordinates,  z  positive  in  the 
upstream  direction  and  p  positive  outward  from  the  axis  of 
symmetry,  E,  and  n  are  oblate  spheroidal  coordinates,  c  is  a 
positive  constant,  and  i  =  .  The  azimuthal  coordinate,  <j> , 

is  the  same  in  both  systems.  Since  c  >  0,  this  leads  to  the 
relations : 

z  =  c  sinhC  cosn  (4. 2. 1.2) 

p  =  c  coshC  sinn  (4. 2. 1.5) 

Each  coordinate  is  limited  by  the  ranges: 

0  <  i  ±  0  £  1  1  ^  >  0<.<)><  2it  (4. 2. 1.4) 

Referring  to  Figure  4. 2. 1.1,  surfaces  of  constant  £  are 
oblate  spheroids,  so  the  geometry  is  suited  for  selecting  c 
such  that  £  is  constant  on  the  surface  of  the  spheroid.  By 
defining  the  axis  ratio  AR  =  b/a,  where  a  is  the  major  axis 
and  b  is  the  minor  axis,  the  following  expressions  hold: 

b  =  c  sinh  5  ,  n  =  0  (4. 2.1. 5) 

a  =  c  cosh£Q,  n=  tt/ 2  (4. 2. 1.6) 

and  at  the  surface  of  constant  it  holds  that 

b/a  =  tanh  50  =  AR  (4.2.17) 

which  leads  to  £0  defined  by: 
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«o  '  7  l0«  (f-Hof)  (4. 2. 1.8) 

7  £ 

Happel  and  Brenner13  review  the  metric  coefficients  and 

give  the  expressions  needed  to  derive  the  various  operators 

of  vector  calculus.  The  most  important  one  for  this  study  is 
2 

the  E  operator  in  oblate  spheroidal  coordinates  given  by: 

e2  .  (sinh2c  ♦  cos2n).  ,iL  ♦  » 2  tanh?  5 

c2  dC  3n2  H 

-  cotn  £-}  (4. 2. 1.9) 


where  c  is  given  by  Equation  4. 2. 1.6. 

4.2.2  Circular  Cylindrical  Coordinates  and  Transformations 

Since  this  coordinates  system  is  widely  used  in  engi¬ 
neering,  it  will  not  be  reviewed  in  detail.  However,  trans¬ 
formation  relationships  from  oblate  spheroidal  to  cylindrical 
coordinates  are  needed.  From  Happel  and  Brenner3**  it  can  be 
shown  that  the  transformation  of  vector  components  from  oblate 
spheroidal  to  cylindrical  coordinates  is  given  by: 

u^  cosh£  cosr)  -  u^  sinh£  sinn 
Uz  (sinhZ£  +  cosZn)1/Z 

u^  sinh£  sinn  +  cosh£  cosn 

UP  =  -  "  uZe - TTU2 - 

(sinh  £  +  cos  n) 


(4. 2. 2. 2) 
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u 


(4. 2. 2. 3) 


where  the  vector  components  are: 

u(p,*,z)  =  upip  +  ♦  uziz  (4. 2. 2.4) 

uU.Kn)  -  Vc  +  +  unin  (4. 2. 2. 5) 

for  the  cylindrical  and  oblate  spheroidal  coordinates, 
respectively. 

Figure  4. 2. 2.1  represents  the  collision  between  an 
oblate  spheroid  and  a  sphere.  A  sphere  is  said  to  have  col¬ 
lided  with  the  oblate  spheroid  when  the  sphere  center  was 
within  one  sphere  radius  of  the  spheroid  surface.  Due  to 
the  geometry  of  the  two  particles,  the  scheme  is  involved 
but  still  only  requires  two  dimensions  since  the  oblate 
spheroid  has  rotational  symmetry  about  its  minor  axis. 

The  cylindrical  coordinates  of  the  curve  representing 
the  intersection  of  the  oblate  spheroidal  surface  with  the 
meridional  plane  is: 

A  A  A  A 

zxiz  +  pi^p  =  cos$iz  +  sin<J>ip  (4, 2. 2. 6) 

where  the  angle  <p  is  the  angle  which  gives  the  shortest  dis¬ 
tance  between  the  center  of  the  sphere  and  the  surface  of 
the  oblate  spheroid.  Equation  4. 2. 2. 6  has  been  nondimension- 
alized  by  dividing  by  the  semi-major  axis  length,  a,  of  the 
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oblate  spheroid.  Referring  to  Figure  4. 2. 2.1,  the  separation 
distance,  S,  is: 

S  =  ( z 2  "  +  ( P 2  *  Pj) 

2  2  2  2  2 
=  ( z 2  +  P2  +  AR  cos  <p  +  sin  <p 

-  2ARz ^cosij)  -  2p2cos<J>  (4. 2. 2. 7) 

Since  S  is  a  function  of  the  angle  <|>,  the  closest  distance 
is  achieved  when  dS/d4»  is  zero.  Differentiating  and  setting 
the  result  equal  to  zero  gives: 

2 

sin<pcosip(l  -  AR  )  +  ARz2sin<J>  -  p2cos<j>  =  0  (4. 2. 2. 8) 

The  angle  0  is  easily  determined  by: 


so  that  the  angle  <p  can  be  determined  implicitly  by  the  fol¬ 
lowing  expressions: 

f ( 4> )  =  sin<t>(l  -  AR2)  -  Z2(tan0  -  ARtan<{>)  =  0  (4.2.2.10) 

0°<  0  <  45° 

f(<$>)  =  cos<f>(l  -  AR2)  -  p2(cot<|>  -  ARtane)  *  0  (4.2.2.11) 
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Both  equations  are  stable  and  converge  rapidly  to  the  angle 
4>  for  a  sufficiently  good  initial  estimate  is  made.  In  an 
iteration  scheme,  the  initial  guess  should  be  <p  =  0  and  then 
allow  the  scheme  to  converge  to  an  error  <  10~14  or  smaller. 

If  the  Newton- Raphson  method  of  finding  the  zeros  of  a 
function  is  used,  convergence  is  very  fast.  Derivatives  are: 

7  Z?AR 

f'  (<|j)  =  co  s  C 1  -  Air)  +  -=—j-  =  0  (4.2.2.12) 

cos  *<p 

0°  <  Q  <  45° 

f "* ( 4> )  =  sin4> (1  -  AR2) - =  0  (4.2.2.13) 

sin  <)> 

45°  <  9  <  90°. 

Again  the  initial  guess  should  be  <J>  =  0. 

4 . 3  Viscous  Flow  Past  Oblate  Spheroids 

The  steady,  axisymmetric  incompressible  flow  around  an 
oblate  spheroid  is  described  by  the  Navier- Stokes  equation. 
This  equation  is  nonlinear,  second-order  partial  differential 
equation  with  three  unknown  quantities,  i.e.,  two  velocity 
components  and  pressure.  Pressure  is  eliminated  by  taking 
the  curl  twice  of  the  Navier-Stokes  equation,  and  the  Stokes 
stream  function,  which  is  related  to  the  velocity,  reduces 
the  unknown  velocity  components  to  one  unknown  quantity. 
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This  results  in  a  fourth-order  partial  differential  equation 
with  only  the  stream  function  unknown.  By  introducing  vortic- 
ity  as  the  curl  of  velocity,  this  fourth-order  differential 
equation  is  split  up  into  two  second-order  coupled  partial 
differential  equations.  This  procedure  is  well  known  and  can 
be  found  in  standard  fluid  mechanics  textbooks. 

Using  oblate  spheroid  coordinates  and  the  above  proce¬ 
dure,  the  Navier- Stokes  equations  can  be  written  as: 

E 2  if/  -  G  (4.3.1) 

7  Re  cosh?  cosh?  sinn 

EZG  =  - t-2 - - -  J  OKF)  (4.3.2) 

2 (sinh  ?  +  cos^n)  ^,n 

where : 

ip  ~  Stokes  stream  functions  ‘ 

G  ~  modified  vorticity  defined  by  Equation  4.3.3 
F  ~  modified  vorticity  defined  by  Equation  4.3.4 
Jj.  ^  ~  Jacobian  operator  defined  by  Equation  4.3.6 
Re  ~  Reynolds  number. 

The  two  modified  vorticities  are  defined  in  terms  of  the 
dimensionless  vorticity  u>  as: 


cosh?sinn 

U1  - TTT - - 

cosh? 

o 

ucosh?0 
cosh?  sinn 


(4.3.3) 


F 


(4.3.4) 
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where  to  is  defined  in  terms  of  the  velocity  ^  by: 

Z  =  curl  $  (4.3.5) 

the  Jacobian  operator  is  simply: 


(*,F) 


dtp  3F  _  dW  3F 

sT  "5TT  ITn  "5iF. 


(4.3.6) 


The  velocity  is  related  to  the  stream  function  and  is  given 
in  oblate  spheroidal  coordinates  as: 


cosh  £ 


uc  -  - 


o  an 

T 


cosh£  sinn  (sinh^E;  +  cos2n)1^ 


(4.3.7) 


cosh^C 


u 


dtp 
o 


2  j  1  /  7 

cosh£  sinn  (sinh  E;  +  cos  n)  ' 


(4.3.8) 


%  =  ufir  +  u  i  .  (4.3.9) 

5  C  n  n  v  ' 

The  following  dimensionless  quantities  were  used: 

<j>  =  ^*/uroa2  (4.3.10) 

co  =  ato*/uoo  (4.3.11) 

2au 

00 


Re 


v 


(4.3.12) 
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where  u®  is  the  velocity  of  the  undisturbed  stream,  v  is  the 
kinematic  viscosity  of  the  fluid,  ij/*  is  the  dimensional 
stream  function  and  to*  the  dimensional  vorticity.  The  quan¬ 
tity  a  is  the  length  of  the  semi-major  axis  of  the  spheroid, 
and  is  equal  to  c  cosh^Q,  where  E,Q  is  given  by  Equation 
4. 2. 1.8. 

Boundary  conditions  are  based  on  the  Stokes  stream  func¬ 
tion  and  vorticity.  As  the  partial  differential  equations 
are  of  second  order,  four  boundary  conditions  are  needed  to 
uniquely  determine  a  solution.  They  are: 

Along  the  axis  of  symmetry:  n  =  0,tt 

=->  \p  =  0 ;  us  =  0  (4.3.13) 

Along  the  surface  of  spheroid:  £  =  E,q 

=s>  ip  =  0 ;  w  =  G/sinn  (4.3.14) 

At  the  outer  boundary:  £  = 

=  i^b ;  ui  =  0  (4.3.15) 

2  7 

,  sin2  cosh“5u 

where  4>u  =  j  - 2 -  and  ^b  *s  an  outer  spheroid  shell 

cosli£n 

where  ur/u  <  0.999  at  p  =  0 . 

The  solution  to  Equations  4.3.1  and  4.3.2  for  the  case 
of  Stokes  flow  past  an  oblate  spheroid  has  been  discussed  by 
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Happel  and  Brenner-3  based  on  the  work  of  Oberbeck-3  .  For 
cases  where  Re  >  0 ,  numerical  methods  must  be  used  and  Chap- 

> 

ter  V  will  discuss  these  solutions. 

4 . 4  Equations  of  Aerosol  Motion 

4.4.1  Derivation  of  Dimensional  and  Nondimensional  Equations 

The  equations  of  aerosol  particle  motion  are  derived 
by  summing  the  external  forces  acting  on  a  particle.  Typical 
forces  include  gravitational,  hydrodynamic,  electrostatic, 
thermphoretic  and  diffusiophoretic .  To  determine  the  gravita¬ 
tional  collision  efficiency  in  this  study,  only  gravitational 
and  laminar  shear  forces  are  considered.  The  equation  of 
particle  motion  is: 

d$. 

<mi>  ^  =  [mil]  ♦  (4. 4. 1.1) 

where : 

<mi>  -  mass  the  it^1  aerosol  particle,  allowing  for 
♦  added  mass  effects 

[n^g]  ~  gravitational  force,  allowing  for  fluid 
buoyancy  effects 

Vi  ~  velocity  of  the  it*1  aerosol  particle 
F^  -  drag  force  of  the  i**1  particle. 

The  added  mass  effect  of  the  fluid  as  the  particle  moves 
through  the  fluid  cannot  always  be  neglected  since  the  added 
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(4. 4. 1.7) 
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+ 


t 


and  the  error  introduced  in  Equation  4. 4. 1.1  if  neglected  is: 


Error 


1 


(1 


(4. 4. 1.8) 


From  Table  4. 4. 1.1  it  is  apparent  that  for  spheres,  both 
effects  can  be  neglected  since  the  error  is  on  the  order  of 
10  .  However,  as  the  AR  +  0,  the  strong  effect  of  added 

mass  dictates  that  this  term  be  retained  in  Equation  4.4. 1.1. 
Thus  there  is  no  justification  for  neglecting  the  added  mass 
term  and  retaining  buoyancy  term,  since  buoyancy  is  only 
negligible  if  the  added  mass  is  negligible. 

Equation  4. 4. 1.1  can  be  expanded  into  six  to  eight 
dimensional  equations  depending  on  whether  relative  or  abso¬ 
lute  position  is  needed.  For  this  research  project,  only 
six  equations  are  used  and  are  written  as: 


<m  >  dVzl  =  [m,g]  +  F  , 
L  dt  1  Zi 


dVp 
<ml>  "3F 


=  F 


pl 


*  <m2>  d^z2  =  [m2g]  +  Fz2 


*  <m7>  dVp2  =  F  7 

2  -at"  p2 


*  =  V  -  V 

dt  Vz  2  vzl 


(4. 4. 1.9) 


(4.4.1.10) 


(4.4.1.11) 


(4.4.1.12) 


(4.4.1.13) 
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Table  4. 4. 1.1 


Added  Mass  and  Buoyancy  Effects 


Added  Mass  Effect 


AR 

aa 

Error 

0.001 

636 

7.3(101) 

0.01 

63.5 

2.1(10'1) 

0.05 

12.5 

5.1(10'2) 

0.10 

6.2 

2.6(10’2) 

0.50 

1.1 

4 . 7 ( 10 ’ 3) 

0.75 

0.70 

3 . 0 (10 ’ 3) 

1.00 

0.50 

2 . 1 (10  *  3) 

Buoyancy  Effect 


Error 


4.3(10*3) 


tTable  Based  on  Following  Densities 
pf  =  1.29(10’3)g/cm3,  pr  =  3 . 0 (10’ :) g/ cm3 
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*  {[§■  =  yP2  *  VPl  (4.4.1.14) 

•  The  asterisk  before  the  equation  means  that  all  quanti¬ 

ties  carry  dimensional  units  whereas  if  the  asterisk  appears 
over  individual  quantities  then  only  those  quantities  carry 
dimensional  units.  Equations  4. 4. 1.9  -  4.4.1.14  assume  that 
the  oblate  spheroid  (particle  one)  has  a  larger  terminal 
velocity  than  the  sphere  (particle  two). 

The  expressions  for  the  mass  of  the  oblate  spheroid  and 
sphere  are  given  by  Equations  4.4.1. 5.  The  semi-major  axis, 
a,  is  related  to  the  volume  equivalent  radius  of  the  spheroid 
by : 


* 


a  = 


(4.4.1.15) 


To  render  the  six  equations  of  motion  nondimensional , 

the  nondimensionalizing  unit  of  length  is  taken  as  the  volume 

equivalent  radius  of  the  oblate  spheroid,  r.  ,  and  the  non- 

f  v 

dimensionalizing  unit  of  velocity  is  taken  as  the  Stokes 
terminal  velocity  of  the  larger  particle,  V^,  where: 


2  plg  rl,v 
9  bf 


(4.4.1.16) 


Again,  fluid  properties  are  viscosity,  and  density,  P£. 

Time  is  nondimensionalized  with  respect  to  particle  one: 
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*  ‘l 

rl,v 

=  rf 

00 1 

(4.4.1.17) 

Each 

of  the  variables  in  the 

dimensional  equations, 

Equations 

4. 4. 1.9  -  4.4.1.14,  is  nondimensionalized  as: 

Vzl 

-  V^/V*! 

(4.4.1.18) 

VP1 

*  vjj/v:. 

(4.4.1.19) 

Vz2 

•  Vh/VJl 

(4.4.1.20) 

V  -7 

p2 

-  vPVv£i 

(4.4.1.21) 

Z 

■  zVrl,v 

(4.4.1.22) 

P 

-  p‘/rf>v 

(4.4.1.23) 

t 

=  t*/tj 

(4.4.1.24) 

The  force 

terms  are  written  in  the 

* 

general  form  as : 

Fn 

■  Fji  ri,vvii 

(4.4.1.25) 

rn 

-  Fj2*/6nuf  rJ^V^ 

(4.4.1.26) 

where  j  = 

z  or  p,  the  cylindrical 

coordinates . 

The  nondimensional  groups  and  four  nondimensional  con 


stants  are  now  defined.  The  first  constant,  r,  is  the  ratio 
of  the  particle  volume  equivalent  radii, 
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r  ■  r5,v/r*,v 


(4.4.1.27) 


the  second,  y,  is  the  ratio  of  the  particle  densities, 


Y  =  P^/pf 


(4.4.1.28) 


the  third,  B,  is  the  ratio  of  the  particle  buoyancy  effects, 


a  -  ;f/Sp 

D  ~  - * - JTT— 


(1  *  Pf/Pj) 


(4.4.1.29) 


and  the  final  constant,  I,  is  the  particle  inertial  effects 
due  to  added  mass, 


*  * 


1  +  (Pf/Pp  aa„ 


1  +  (Pf/Sp  A, 


(4.4.1.30) 


so  that  the  first  nondimensional  group,  L,  relates  the  dif¬ 
ferences  between  the  particle  radii,  density,  and  shape,  that 


is : 


L  = 


r4y2  B I 


r2  v  4  p2  2 

-  (4^)  (4)  c 


*  .  * 


*  * 


i  -  pf/p 

( 


l  +  (pf/pp  aa 


l,v 

.  ,iL* 


i  -  pf/p{ 


i  +  (Sf/Sp  V 


-) 


( 

r 


* 

2  -  pf 

)  (4 - 4)  (tt 


* 

4  p: 


Po  +  Pf  A 


f  A. 


l.v 


P1  '  pf 


pl  +  pf  AA, 


(4.4.1.31) 


I 

1 


69 


The  second  nondimensional  group  is  the  Stokes  number,  Stk, 
and  is  the  ratio  of  the  penetration  distance  of  the  small 
particle  to  the  characteristic  radius  of  the  obstacle  (pro¬ 
jected  radius  of  oblate  spheroid  for  the  case  being  con¬ 
sidered).  If  the  smaller  particle  were  injected  into  motion¬ 
less  fluid  with  a  velocity  equal  to  its  Stokes  terminal  veloc¬ 
ity,  Vot2»  the  particle  will  come  to  rest  after  travelling  z2, 
its  penetration  distance.  Integration  of  the  equations  of 
motion  gives  the  penetration  distance  as: 


* 


Z2 


<m~> 

- 1 -  v 

^fr2,v  ~2 


and  the  Stokes  number  is: 


(4.4.1.32) 


Stk  =  z*/a*  =  AR1/3z*/r1*y. 
The  Stokes  number  is  also  given  by: 


* 

V 


(1  +  (Pf/Pp  Aa  ) 

Stk  -  -i±*  r\2B  - r-, - 

g  a  (1  -  Pf/  Pi') 


(4.4.1.33) 


Thus  the  six  nondimensional  equations.  Equation  4. 4. 1.9  - 
4.4.1.14,  are  now  written  as: 


dVzl 

-ar-  ■  AR 


1/3  stir  o  *  Fzi> 


(4.4.1.34) 


gjUl 
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dV 


Pi  _  AD 1 / 3  L  D 

ar-  _  AR  stir  Fpi 


dVz2  _  AR1/3  L  fR  Fz2. 

T  Stic  (B  '7?) 


“3T 


dVp2 

~cTt~ 


AR1/3  L  rFp2s 

T  Stic 


yr 


dz 

3t 


Vz2  -  Vzl 


=  v  -  v 

<Tt  Vp2  vpl 


4.4.2  Drag  Force  Terms:  Superposition  Method 


(4.4.1.35) 

(4.4.1.36) 

(4.4.1.37) 

(4.4.1.38) 

(4.4.1.39) 


The  drag  force  in  Equations  4.4.1.34  -  4.4.1.39  must  be 
known  before  these  equations  can  be  integrated.  As  previously 
mentioned,  no  analytical  expressions  are  available  to  model 
the  hydrodynamic  interaction  between  an  oblate  spheroid  and 
a  sphere.  It  is  therefore  necessary  to  use  the  so-called 
superposition  scheme.  In  this  method,  the  motion  of  each 
particle  is  superpositioned  onto  a  fluid  velocity  field.  For 
the  superposition  method,  the  nondimensional  force  terms  are: 


CnRe 
<-§?->  i 

* 

*  * 

(4.4. 2.1) 

jl  = 

(vji 

■  V/v-i 

CnRe 

<- vrh 

* 

*  * 

(4. 4. 2. 2) 

n  = 

CVj2 
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where  j  s  z,p  and  CD  is  the  drag  coefficient  for  the  itn  par¬ 
ticle.  The  velocities  of  the  particles,  V j ^  and  Vj2j  are 
obtained  from  integration  of  the  equations  of  motion.  The 
fluid  velocity  terms  caused  by  motion  of  the  two  particles 
are  W-.  and  W.~.  In  order  to  evaluate  the  nondimensional 

J  i  J4 

force  terms,  the  flow  fields  around  the  oblate  spheroid  and 
sphere  must  be  known.  Several  possibilities  exist  which 
include  using  Stokes  solution,  Oseen  solution,  or  a  numerical 
approximation  to  the  full  Navier-Stokes  equation. 

Since  the  oblate  spheroids  were  assumed  to  be  the  col¬ 
lecting  particle,  it  was  decided  to  solve  the  full  steady- 
state  Navier-Stokes  equation  as  described  on  page  58  .  The 
possibility  of  having  Reynolds  numbers  greater  than  one  is 
quite  likely.  On  the  other  hand  for  very  small  spheres,  the 
particles  being  collected,  there  is  some  advantage  in  using 
an  analytical  expression  for  the  velocity  field.  However, 
if  the  Reynolds  number  is  large  enough  to  invalidate  Stokes 
solution,  the  velocity  field  can  be  calculated  by  using  an 
oblate  spheroid  with  an  axis  ratio  0.999.  For  Stokes  flow, 
the  analytical  expression  for  the  velocity  field  around  a 
sphere  in  cylindrical  coordinate  system  is: 


V  -  c<* 


r2 , v  pz 
(PZ  +  zZ) 


a  — riyT)} 

(pZ  +  zZ) 


(4.4.2.3) 


These  expressions  are  valid  for  Re  <<  1.  For  larger  Reynolds 
numbers,  it  is  appropriate  to  solve  the  full  steady- state 
Navier-Stokes  equation. 


CHAPTER  V 


NUMERICAL  METHODS  AND  COMPUTER  CODE  DEVELOPMENT 


5.0  Introduction 

In  the  previous  chapters,  the  Navier-Stokes  equation, 
the  equations  of  aerosol  particle  motion,  and  other  relation¬ 
ships  are  derived  in  order  to  determine  the  gravitational 
collision  efficiency.  This  chapter  explains  the  numerical 
methods  used  and  analyzes  the  problems  associated  with  each. 
These  methods  are  used  in  the  computer  code  NGCEFF,  a  code 
specifically  tailored  to  calculate  the  gravitational  colli¬ 
sion  shape  factor,  3,  for  oblate  spheroidal  shaped  particles. 

The  last  portion  of  this  chapter  explains  how  the  code 
works  and  various  options  available.  The  code  is  written  as 
a  series  of  subroutines  and  function  routines,  thereby  allow¬ 
ing  for  modifications  and  improvements.  This  is  done  so  that 
other  nonspherical  shaped  particles  can  be  substituted  with¬ 
out  writing  a  new  code  from  the  beginning.  Hopefully  the 
comment  statements  and  descriptive  variable  names  will  aid  in 
this  effort. 

5 . 1  Numerical  Methods  and  Analysis 


5.1.1  Solution  of  the  Equations  of  Motion 
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5. 1.1.1  Numerical  Integration 

The  solution  of  the  six  nondimensional  equations,  Equa¬ 
tions  4.4.1.34  -  4.4.1.39,  require  that  an  accurate  numerical 
scheme  be  used.  The  problems  involved  in  integrating  a  set 

of  equations  exhibiting  the  mathematical  property  of  stiff- 

38  12 

ness  are  discussed  by  Gear  and  also  Pertmer  ,  who  addressed 

this  question  specifically  for  the  problem  at  hand. 

A  stiff  set  of  differential  equations  can  be  defined  as 

a  set  of  differential  equations  which  have  eigenvalues  with 

negative  real  parts  that  are  widely  separated  from  each  other. 
38 

From  Gear  ,  consider  the  following  set  of  equations: 


u'  *  998u  +  1998v 

v'  =  - 999u  -  1999V  (5. 1.1. 1.1) 


This  system  has  two  eigenvalues,  -1  and  -1000,  and  its  solu¬ 
tion  is: 


u 

v 


,  -t  -lOOOt 

2e  -  e 

- 1  -lOOOt 

-e  +  e 


(5. 1.1. 1.2) 


Although  the  terms  e'1000t  die  away  rapidly  with  increas¬ 
ing  time  t,  they  dominate  the  solution  initially  and  control 
the  stability  and  accuracy  of  the  numerical  method.  Decreas¬ 
ing  the  step  size  during  integration  will  help  if  enough  sig¬ 
nificant  digits  are  available.  This  significantly  increases 
the  time  of  integration,  and  if  the  integration  routine  is 
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part  of  an  inner  scheme  which  in  turn  is  part  of  an  outer 
iteration  scheme,  then  this  solution  may  be  impractical.  If 
not  enough  significant  digits  are  available  for  the  machine 
being  used,  then  any  solution  is  really  meaningless  since 
excessive  errors  will  be  introduced  by  truncation  or  roundoff 

procedures  in  the  initial  steps. 

38 

Gear  in  his  book  defines  stiff  stability  and  proceeds 

to  determine  the  existence  of  stable  systems  of  order  up  to 

six.  Basically,  Gear's  method  is  a  modification  of  an  Adams 

39 

multistep  predictor-corrector  method.  The  IMSL  subroutine, 
DGEAR,  makes  use  of  Gear's  subroutine  DIFSUB.  A  complete 
discussion  of  Gear's  method  can  be  found  in  his  book,  and 
DGEAR  is  available  from  IMSL  or  many  computing  centers  at 
universities . 

12 

From  the  numerical  calculations  of  Pertmer  ,  the  prob¬ 
lem  of  mathematical  stiffness  gets  increasingly  difficult  as 
the  particle  sizes  decrease.  His  investigation  of  the  vari¬ 
ables  in  the  aerosol  equations  of  motion  shows  that  the  most 
rapidly  changing  variable  was  a  particle's  horizontal  veloc¬ 
ity.  This  is  a  reasonable  result  and  certainly  one  which 
applies  to  the  present  case.  Pertmer  believes  that  the 
stiffness  of  the  aerosol  equations  of  motion  is  due  to  large 
horizontal  accelerations  of  both  particles. 

Since  Gear's  method  successfully  solved  the  six  non- 
dimensional.  equations  of  motion  derived  by  Pertmer,  it  was 


«  • 


natural  to  employ  this  method  for  the  six  equations  derived 
in  Chapter  V.  The  initial  conditions  needed  to  integrate 
these  equations  are  based  on  the  Stokes  terminal  velocity  and 
the  volume  equivalent  radius  of  particle  one.  It  is  assumed 
that  the  aerosol  particles  are  falling  vertically  downward 
with  velocities  equal  to  their  actual  terminal  velocities 
based  on  their  Reynolds  numbers.  The  terminal  velocities 


9  p  B  1,V 


(1  *  P/Pi") 


*  V®  =  Re  y/2ap 

cl  a 


*  V“v,l  '  Rev.l  p/2rl,vp 


V“v,2  '  I  v  8  r22,v  n  •  p/pP 


(S. 1.1. 1.3) 


(5. 1.1. 1.4) 


(5. 1.1. 1.5) 


(5. 1.1. 1.6) 


where : 


~  Stokes  terminal  velocity  of  particle  one  based 
on  its  volume  equivalent  radius,  r. 

1  ,  V 

~  Terminal  velocity  of  particle  one  based  on  its 
semi-major  axis  length,  a 
~  Terminal  velocity  of  particle  one  based  on  its 
volume  equivalent  radius,  r^  y 
~  Stokes  terminal  velocity  of  particle  two  based 
on  its  volume  equivalent  radius,  ^  v* 
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In  Chapter  IV,  the  equations  of  motion  were  nondimension- 
alized  with  respect  to  and  r^  y.  Assuming  the  initial 
vertical  separation  of  the  particles  is  taken  to  be  fifty 
times  the  volume  equivalent  radius  of  the  larger  aerosol  par¬ 
ticle,  the  initial  horizontal  velocities  are  set  equal  to 
zero.  The  initial  conditions  now  become: 


V 


zl) 


o 


V  °°a/V  °°1 ;  v  %  j/V  00 1 


(5. 1.1. 1.7) 


V 


z2) 


o 


(5. 1.1. 1.8) 


V,0  •  ° 


V  „  =  0 

P2)0 


z.  =50 
'  o 


(5. 1.1. 1.9) 


(5.1.1.1.10) 


(5.1.1.1.11) 


The  initial  vertical  separation  is  somewhat  arbitrary.  It 
should  be  large  enough  so  that  one  aerosol  particle  has  little 
or  no  effect  on  the  initial  motion  of  the  other  aerosol  par¬ 
ticle  but  small  enough  so  that  integration  time  is  appropriate 

1 2 

and  worthwhile.  Pertmer  found  that  50  was  a  satisfactory 

35 

compromise,  as  did  Pitter  and  Pruppacher  in  their  investi¬ 
gation.  » 
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DEC  80  R  F  TUTTLE 

UNCLASSIFIED  AFIT-C 1-80-820  _  _  _ 
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The  aerosol  particle  initial  horizontal  separation,  p.  , 

•'o 

is  the  parameter  which  determines  the  gravitational  collision 
efficiency.  From  the  integration  of  the  six  equations  of 
aerosol  motion,  the  critical  initial  horizontal  separations 
between  the  two  aerosol  particles  can  be  determined.  From 
the  results  of  Pitter  and  Pruppacher^5 ,  the  gravitational 
collision  efficiency  must  consider  the  possibility  that  the 
collision  domain  is  a  circular  annulus.  Thus  the  definition 
needs  to  be  expanded  to: 

(5.1.1.1.12) 

where  Epc  is  the  sum  of  collision  domains  above  the  particle 

35 

in  question.  Pitter  and  Pruppacher  found  only  one  colli¬ 
sion  domain,  which  is  defined  as  the  difference  between  the 
largest  and  smallest  initial  horizontal  separations  which 
give  grazing  collisions.  Collisions  were  always  assumed  to 
occur  if  the  initial  horizontal  separation  was  greater  than 
the  smallest  and  less  than  the  largest  critical  grazing 
values . 

The  possibility  of  having  an  annular  collision  domain 
presents  a  problem  in  any  search  routine.  The  bisection 
method  is  not  suitable  since  cases  can  be  formulated  in  which 
the  scheme  will  fail.  Normally  detail  information  is  not 
available,  so  that  any  elaborate  bisection  scheme  creates 


e  =  (F - - ) 

rl,v  r2,v 
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needless  complications.  The  most  prudent  scheme  is  to  locate 
the  maximum  p  for  a  grazing  collision  and  then  decrease  p  in 
hopes  of  finding  any  "breaks"  or  discontinuities  between  the 
maximum  p  and  p  equal  to  zero. 

5.1.1. 2  Determination  of  Minimum  Separation 

After  each  step  of  integration,  the  oblate  spheroid  moves 
closer  to  the  spherical  particle.  Associated  with  each  step 
is  an  angle  of  separation,  <j>,  and  the  corresponding  separation 
distance,  S  (see  Figure  4. 2. 2.1).  If  a  direct  or  grazing 
collision  occurs,  the  integration  routine  will  stop.  However, 
if  the  smaller  particle  misses  the  larger,  nonsperhical  par¬ 
ticle,  the  minimum  separation  does  not  have  to  correspond  to 
one  of  the  separation  distances  calculated  for  each  time  step. 
Figure  5. 1.1. 3  illustrates  what  may  occur  during  integration. 

To  estimate  the  minimum  separation  distance,  S  .  ,  the  mini- 
mum  separation  angle,  4>ra^n»  is  determined  by  fitting  a 
natural  cubic  spline  to  the  data  points.  From  the  cubic 
spline  function,  the  angle  of  minimum  particle  separation  is 
found  by  taking  the  derivative  of  the  spline  and  setting  it 
equal  to  zero.  Of  the  two  solutions,  only  one  angle  will 
lie  between  <J>j  and  <p To  find  minimum  particle  separation, 

the  natural  cubic  spline  calculates  it  by  using  Spline 

39 

coefficients  are  determined  by  using  IMSL  subroutine  ICSICU 
and  interpolation  is  accomplished  using  IMSL  subroutine  ICSEVU. 
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5.1.2  Numerical  Methods  Used  to  Calculate  Velocity  Fields 

The  equations  of  aerosol  motion  require  that  the  veloc¬ 
ity  fields  around  both  particles  be  known  in  order  to  esti¬ 
mate  the  drag  forces  (see  Equations  4. 4. 2.1  and  4. 4. 2. 2). 

The  flow  field  around  the  smaller  particle  is  obtained  from 
Stokes  analytical  expression.  The  flow  field  around  the 
larger  particle,  i.e.,  either  the  oblate  spheroid  or  its 
volume  equivalent  sphere,  must  be  calculated  from  the  solu¬ 
tion  of  the  Navier-Stokes  equation.  The  solution  of  this 
equation  is  in  terms  of  the  stream  function  and  vorticity. 

To  find  the  flow  fields,  Stokes  relationships  must  be  used 
with  some  appropriate  numerical  method.  By  knowing  the  loca¬ 
tion  of  the  center  of  each  interacting  particle,  the  velocity 
field  of  each  can  be  interpolated  at  the  location  of  the 
other  particle.  Thus  two  numerical  procedures  must  be  used. 
The  first  has  to  do  with  the  solution  of  Stokes  relationships 
and  the  second  is  the  interpolation  of  the  velocity  fields  at 
the  proper  locations. 

5. 1.2.1  Numerical  Solution  of  Stokes  Relationships 

From  Equation  4.2. 2.1  and  4. 2. 2. 2,  the  velocity  field  in 
cylindrical  coordinates  is  easily  obtained.  Determination  of 
oblate  spheroidal  velocities  and  is  more  difficult. 
Formally  and  are  given  by  the  following  expressions: 
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cosh2  £  3^/3n 

U  - ^-7 - 2 T7T  (5. 1.2. 1.1) 

^  coshc  sinn  (sinh  £  +  cos  n) 

cosh2  £  3i^/3C 

U  =  - -2-1 - 2__  (5. 1.2. 1.2) 

n  cosh£  sinn  (sinli  £  +  cos  n) 


4 


I. 

1 

I 

I 


27  35 

Previous  investigators  *  used  a  finite  difference  scheme 
with  second  order  accuracy  to  compute  the  velocity  at  each 
grid  point. 

It  was  decided  that  piecewise  cubic  splines  offered  cer¬ 
tain  advantages  over  the  finite  difference  scheme.  First, 
since  cubic  splines  are  derived  by  matching  the  first  and 
second  derivatives  (and  thus  the  slope  and  curvature)  at  each 
grid  point  in  order  to  form  a  piecewise  continuous  (smooth) 
curve,  the  fitting  of  the  stream  function  field  with  cubic 
splines  gives  directly  the  derivatives  at  each  grid  point. 

All  that  is  required  to  calulate  the  partial  derivatives  in 
Equations  5. 1.2. 1.1  and  5. 1.2. 1.2  is  to  take  the  derivative 

of  a  spline,  which  is  simply  the  derivative  of  a  polynomial 

40 

of  order  4  or  of  degree  3  (deBoor  ) . 

Since  two  partial  derivatives  are  required,  two  sets  of 
splines  are  formed.  One  set  contains  those  splines  con¬ 
structed  with  respect  to  the  coordinate  n  (holding  £  constant) 
and  the  other  set  contains  those  splines  constructed  with 
respect  to  the  coordinate  £  (holding  n  constant).  Splines 
and  their  corresponding  spline  coefficients  were  computed 


diiikMMiaibA 


Til  ''ifn 
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to 

using  the  IMSL  routine  ICSICU.  ICSICU  computes  the  piece- 
wise  continuous  splines  coefficients  when  given  a  set  of 
points  and  arbitrary,  but  somewhat  appropriate,  second  deri¬ 
vative  end  conditions.  From  the  IMSL  description  of  ICSICU, 
define  f^",  ^2"’  ^k- 1 '  anc*  ^k"  w^ere  ^  t^ie  number  of  data 
points,  as  the  second  derivatives  of  the  curve  formed  by  the 
set  of  data  points.  The  arbitrary  second  derivative  end  con 
ditions  are: 


+  <hf2"=  q2 


(5. 1.2. 1.3) 


3  k-1 


+  2  f, 


q4 


(5. 1.2. 1.4) 


where  q^,  i  =  1,...,4  are  the  four  arbitrary  input  boundary 
conditions . 

Equations  5. 1.2. 1.3  and  5.1. 2.1.4  point  out  the  second 
advantage  of  splines  over  the  finite  difference  method.  Com¬ 
putation  of  on  the  axis  of  symmetry  presented  Pitter  and 
Pruppacher^5  with  a  problem  because  when  n  =  0 ,  both  the  sinq 
and  3iji/8n  are  zero,  necessitating  the  use  of  L'Hospital’s 
rule.  Pitter  and  Pruppacher  found  that  the  use  of  L'Hospital's 
rule  gave  unsatisfactory  numerical  results.  They  had  to  use 
the  adjacent  derivatives  to  estiamte  U^  at  the  grid  locations 
that  gave  the  indeterminate  results.  Spline,  however,  makes 
use  of  the  fact  that  3^/3n  =  0  at  n  =  0  and  it  by  forcing  the 
spline  at  the  end  points  to  have  first  derivatives  equal  to 
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zero.  This  is  accomplished  by  setting  q^.  =  q  *  1  and: 

q2  =  (1»2  '  ^ i ^  /  C An)  2  (5. 1.2. 1.5) 

q4  =  (4»k  -  ^.^/(An)2  (5. 1.2. 1.6) 

This  results  in  "shaping"  the  splines  at  the  end  points  so 
that  U £  can  be  determined  accurately  at  the  adjacent  grid 
points.  To  compute  at  n  =  0,  L 'Hospital’s  rule  can  now 
be  applied  to  Equation  5. 1.2. 1.1  and  the  end  splines  can  be 
differentiated  twice  to  give  3  \p/B  n.  This  will  give  good 
results,  but,  surprisingly,  better  results  can  be  obtained 
by  using  Gregory-Newton  forward  and  backward  extrapolation 
formula,  i . e . : 


n  =  3(U  -  U  )  +  U 

^n-0  ^2  ^3  **4 

(5. 1.2. 1.7) 

n  =  3(U  -  u  )  +  U 

^k  Vl  ^k-2 

(5. 1.2. 1.8) 

Calculation  of  requires  a  different  approach  to 
estimate  the  end  point  parameters.  First  and  second  deriva¬ 
tives  are  not  necessarily  known  but  because  of  the  no-slip 
condition  imposed  at  the  surface  of  the  oblate  spheroid,  the 
value  of  is  already  known,  where  £0  is  the  value  of 
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£  at  the  surface  of  the  spheroid.  It  is,  therefore,  not 
important  to  calculate  U  ..  but  U  .  is  needed,  where 

n,S-€0  n)C-5n 

£n  is  the  value  of  £  at  the  chosen  outer  boundary  envelope. 

By  specifying  the  second  derivatives  at  the  end  points,  the 
correct  curvature  of  adjacent  splines  can  be  insured.  This 
is  done  by  setting: 

qx  =  q3  =  0  (5. 1.2. 1.9) 

q2  =  2(2^  -  5^  +  4^  -  *4)/(A£)2  (5.1.2.1.10) 

q4  =  2 (-*k_ 3  +  4*^_2  "  i  +  2*j^)/(A^)  (5.1.2.1.11) 

where  second  order  forward  and  backward  difference  expres¬ 
sions  are  used  to  estimate  second  derivatives.  Equations 
5.1.2.1.10  and  5.1.2.1.11  assume  that  the  angle  n  is  held 
constant.  The  opposite  was  the  case  for  Equations  5. 1.2. 1.5  - 
5. 1.2. 1.8. 

The  accuracy  of  the  methods  discussed  here  could  be 
checked  since  an  analytical  solution  is  available  for 
Reynolds  number  equal  to  zero.  Thus,  Oberbeck's  solution 
was  used  to  verify  the  improvement  offered  by  Equations 
5. 1.2. 1.5  -  5.1.2.1.11  over  that  offered  by  a  finite  differ¬ 
ence  scheme. 
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5. 1.2. 2  Interpolation  of  the  Velocity  Fields 

In  order  to  use  the  velocity  components  determined  by 
Equations  4. 2. 2.1  and  4. 2. 2. 2,  the  location  of  the  sphere 
center  in  oblate  spheroidal  coordinates  must  be  calculated. 
This  information  gives  the  grid  location  with  the  proper 
indices,  which  are  required  in  order  to  select  the  correct 
velocity  components  from  the  two  velocity  field  arrays. 

To  determine  the  angle  p,  the  following  relationship 
is  used: 


2  2 
P  z  _ 

~7 - - 2 —  - 2 - i 

(1  -  cos  n)  cosS 


(5. 1.2. 2.1) 


where  z  and  p  are  the  known  cylindrical  coordinates,  p  is  the 
unknown  oblate  spheroidal  angle,  and  S  is  semi-major  axis 
length,  a,  divided  by  volume  equivalent  radius  of  the  spher¬ 
oid,  r,  v,  and  the  quantity  cosh£Q.  The  cosh^0  is  the 
oblate  characteristic  length  specified  by  Equation  4. 2. 1.6. 
Solving  for  cos  n,  the  following  expression  is  derived: 

-  D2  +  (P2  D2  +  4  S2  z2)1/2 

cos^n  =  - 5 -  (5. 1.2. 2. 2) 

2  S* 


where : 


D2  =  S2  -  (z2  +  p2) 


(S. 1.2. 2. 3) 
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If  the  positive  sign  in  Equation  5. 1.2. 2. 2  is  selected,  the 
angle  n  is: 


n  =  cos’  *  (x) ,  z  ^  0 
n  =  it  -  cos '  ^  (x)  ,  z  <  0 

where : 

x  =  (cos2n)1,/2,  x  £  1 

otherwise : 

x  =  1. 

To  calculate  £  coordinate,  remember  that: 


(5. 1.2. 2. 4) 

(5. 1.2. 2. 5) 


(5. 1.2. 2. 6) 


(5. 1.2. 2. 7) 


z  =  S  sinh£  cosn  (5.1. 2. 2.8) 

so  that  i  is  given  by: 

C  =  sinh  (g  Cosn-*  (5. 1.2. 2. 9) 

where  S  =  a/r^  v  cosh4Q. 

To  determine  the  indices  needed  to  select  the  correct 
velocity  components,  define: 

£  -  (j  -  1)A?  +  C0,  j  =  1,2,...,  n  (5.1.2.2.10) 


n  =  (i  -  l)An 


» 


i  “  in 


(5.1.2.2.11) 
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where  AC  and  An  are  the  step  sizes  for  the  radial  coordinate 
C  and  angular  coordinate  n,  and  n  and  m  are  the  number  of 


steps 


From  Equation  5. 1.2. 2. 9,  the  index  association  with 


radial  coordinate  C  is 


J  =  (C  -  /AC  +  1 


(5.1.2.2.12) 


and  likewise,  for  angular  coordinate  n,  the  index  is 


I  =  n/An  +  1 


(5.1.2.2.13) 


where : 


j  _<  J  <  j  +1 


i  <  I  <  i  +  1 


(5.1.2.2.14) 


(5.1.2.2.15) 


To  determine  the  value  of  the  velocity  component  at 

location  (I,J),  the  four  velocity  components  around  this 

location  are  properly  weighted.  From  the  discussion  about 

Equation  4. 4. 2.1,  the  velocity  at  location  (I,J)  is  W.?, 

J  z 

j  =  z , p .  Thus  at  ( I ,  J)  : 


Wz2  =  wlUz  +  w2Uz  +  W3UZ 

i  ,3  Z  Zi,3+1  3  Z i  +  i >  3 


+  W/iU, 
4  z 


i+ltj+1 


(5.1.2.2.16) 
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W  -  =  w, U  +  w,U  +  w,U 


'p2  ”1  p 


^  Vi  ^  u 

i,j  2  Pi,j+J 


’3“P; 


i  +  l,j 


+  w.U 

4  Pi  +  1 » j  +1 


(5.1.2.2.17) 


where  w.  .  ,  .  are  the  weighting  factors  defined  by  the 

1  >  1“ A  9  •  •  .  ,4 

following  expressions: 

w1=[i-I+l][j-J+l]  (5.1.2.2.18) 

w _  =  [i  -  I  +  1]  [J  -  j  ]  (5.1.2.2.19) 

w_  =  [I  -  i]  [j  -  J  +  1]  (5.1.2.2.20) 

w4  =  [I  -  i][J  '  j].  (5.1.2.2.21) 


The  results  from  Equations  5.1.2.2.16  and  5.1.2.2.17  are  now 
used  in  Equation  4. 4. 2.1. 

5.1.3  Numerical  Methods  of  Solution  of  the  Navier- Stokes 
Equation 


As  a  first  step  to  compute  the  collision  efficiencies 
of  oblate  spheroids  with  spheres,  the  flow  fields  around 
oblate  spheroids  with  different  Reynolds  numbers  and  axis 
ratios  are  needed.  To  determine  the  collision  shape  factor, 
6,  the  flow  fields  around  volume  equivalent  spheres  with 
different  Reynolds  numbers  are  also  needed.  Since  the 
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Navier-Stokes  equation  is  a  second-order,  nonlinear  partial 
differential  equation,  a  numerical  method  of  solution  is 
necessary. 

Chapter  III  reviews  the  literature  which  reported 
numerical  studies  of  the  flow  past  oblate  spheroids  at  low 
and  intermediate  Reynolds  numbers.  These  studies  show  that 
the  finite  difference  method  can  work  but  that  slow  conver¬ 
gence  and  stability  of  the  iteration  method  are  problems 
that  needed  to  be  analyzed. 

Recalling  Equations  4.3.1  and  4.3.2  central  difference 

2 

expressions  can  be  used  on  the  E  and  Jacobian  operators  at 
grid  locations  where  the  Taylor  series  expansion  is  valid. 

If  the  difference  equations  are  solved  for  the  value  of  the 
stream  function  ip  and  vorticity  G  at  the  points  of  expansion, 
the  results  are: 


♦ij1  "  t  161 


*  e4  -  f  GiJ> 


(5. 1.3.1) 


i  ■*  2|3 9  •  •  #  }  m **  1 
j  ~  2 , 3 , .  ,  *  |  n-1 


where : 

ho  -  2 [ (AC) 2  +  (An) 2] / CAC) 2 (An) 2  (5.1. 3. 2) 


ex  =  (1  -  \  n  tanhO/CAO2 


(5. 1.3. 3) 


•  • 


i 


e2  =  (1  +  j  An  tanh£) /  (A£) 
e3  =  (1  -  jAn  cotn)/(An)2 


e4  =  (1  +  j  An  cotn)/(An)' 


f  =  (sinh2£  +  cos2n) /cosh2£ 


(5.1. 3.4) 


(5.1. 3. 5) 


(5. 1.3. 6) 


(5. 1.3. 7) 


-  radial  coordinate  in  oblate  spheroidal  coordinates 
~  angular  coordinate  in  oblate  spheroidal  coordinates 
~  radial  coordinate  step  size 

-  angular  coordinate  step  size 

~  index  associated  with  radial  coordinate,  £,  see 
Equation  5.1.2.2.10 

~  index  associated  with  angular  coordinate,  q ,  see 
Equation  5.1.2.2.11 

~  value  of  i  that  corresponds  with  the  surface  of 
the  oblate  spheroid 
~  iteration  number. 


The  coefficients  e^y  e2,  e3,  e^,  and  f  are  evaluated  at  the 
grid  location  (i,j)  being  considered.  To  determine  the  pro¬ 
per  value  of  £  and  n  ,  Equations  5.1.2.2.10  and  5.1.2.2.11 
are  used  to  generate  the  vectors  and  arrays. 
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The  vorticity  transport  equation,  Equation  4.3.2,  is 
more  difficult  to  expand.  The  problem  is  caused  by  trying 
to  use  central  difference  expansions  on  the  Jacobian  around 
a  grid  point  next  to  the  boundary  where  n  =  0,tt.  For  example, 
if  one  expands  the  Jacobian  using  a  central  difference  expres¬ 
sion,  one  of  the  factors  which  must  be  evaluated  is: 


F 


k 

i+l,j 


Fi 


1 » j 


(5.1. 3.8) 


where  F  is  defined  by  Equation  4.3.4  and  is  related  to  G  by: 


cosh  £  G 

F  =  - - - (5.1. 3.9) 

cosh  (  sin  n 

When  Expression  5. 1.3. 8  is  evaluated  at  the  adjacent  points 
to  n  =  0 ,  tt  (i.e.,  i  =  2,m-l),  the  sin  is  zero  as  well  as  the 
vorticity.  Application  of  L'Hospital  rule  gives  poor 
results.  This  problem  is  avoided  if  a  forward  difference 
expansion  is  used  when  i=2  and  a  backward  difference  expan¬ 
sion  is  used  when  i=m-l.  For  grid  points  with  3  £  i  j<  m-2, 
a  central  difference  expansion  can  be  employed.  This  pro¬ 
cedure  results  in  three  iteration  formulas: 


G2 J  =  K1{dlG2, j -1  +  d2G2 , j  +1  +  d3G3,j 


*  d4Gi.j> 


j  =  2,3,..., 


n-1 


(5. 1.3. 9) 
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t1 


where  G.  .  =  0  (boundary  condition) , 

l ,  j 


gm  ■  i‘vij-1  *  d6Gi,j.i  *  vtij 


*  Vi.ij1 


(5.1.3.10) 


i  =  3,4,...,  m-2 
j  —  2,3,...,  n-1 


Gm- 1 , j  =  E2{d9Grn-l,j-l  +  d10Gm-l , j  +1  +  dllGm-2,j 


+  di?G™  ^  4 

J.  L  171  “  O  y  J 


(5.1.3.11) 


j  2,3,...,  n-1 


where  Gm  j  =  0  (boundary  condition) .  The  coefficients  in 
Equations  5. 1.3. 9  -  5.1.3.11  are: 


Re  •  3iK  . 

hl  '  <ho  - 7TZ  -’i  > 

cosh  £  sin  n 


(5.1.3.12) 


h2  '  {ho 


Re  •  3^^  . 

+  - , - ^4-} 


cosh  g  sin  n 


(5.1.3.13) 


Re  •  <4  ■ 

l,  =  (e?  +  - - *-} 

1  cosh^(£  +  A£)sin  n 


(5.1.3.14) 


. . 
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dll  =  {e4 


Re  •  4(ii>  ,  .  ,  - 

+ _ _ jl-aj r* _ > 

coshZC  sinZ(n  -  An) 


(5.1.3.24) 


Re  •  (ii^  .  ,  -  i  •  n ) 

i  -  { - y-ri — isiitiiL} 

cosh^C  sin^(n  -  2An) 


(5.1.3.25) 


and  Re  is  the  modified  Reynolds  number  defined  by: 


Re  = 


Re  coshg  sinn 
8  (AC)  (An)  secTir 


(5.1.3.26) 


Coefficients  hQ,  e^,  e^,  and  e^  were  defined  earlier. 


5. 1.3.1  Method  of  Solution  and  Its  Analysis 
41 

Varga  discusses  the  basic  iterative  methods  available 
to  solve  matrix  equations  such  as  Equations  5. 1.3.1  and 
5. 1.3. 9  -  5.1.3.11.  Methods  include  the  point  Jacobi,  Gauss- 
Seidel,  and  successive  overrelaxation  and  underrelaxation 
(SOR)  algorithms.  The  SOR  algorithm  offers  all  the  advan¬ 
tages  of  the  Gauss-Seidel  method  plus  the  advantage  of 
accelerating  convergence. 

Woo^  investigated  the  computational  efficiency  of  the 
SOR,  alternating  direction  implicit  (ADI) ,  and  a  version  of 
the  SOR  method  called  the  dominant  eigenvalue  algorithm 
(see  Orback  and  Crowe**3).  He  applied  these  methods  to  the 
problem  of  viscous  flow  around  spheres,  a  problem  similar  to 
the  one  being  considered  here.  He  found  that  the  dominant 
eigenvalue  algorithm  (DEM)  to  be  3-4  times  faster  than  SOR 


7 

and  ADI.  Pitter  et  al .  ,  applied  DEM  to  viscous  flow  around 

oblate  spheroid  and  claimed  that  it  helped  to  accelerate  con¬ 
vergence. 

The  DEM  algorithm  is  based  upon  the  observation  that 
most  iterations  eventually  approach  a  geometric  progression. 
The  iteration  procedure  continues  until  the  largest  eigen¬ 
value  of  iteration  forcing  matrix  dominates  the  solution. 

The  dominant  eigenvalue  can  then  be  used  to  extrapolate  toward 
the  solution  before  resuming  the  iteration  scheme. 

In  practice,  the  dominate  eigenvalue  is  not  calculated. 
The  DEM  algorithm  is  applied  as  follows.  Divide  the  itera¬ 
tion  of  the  loop  parameter  k  into  two  half  iterations  for  the 
index  j,  the  index  associated  with  the  radial  component  £. 

The  scheme  is  to  use  successive  substitution  with  eigenvalue 
promotion  for  half  of  the  grid  points  and  no  update  for  the 
other  half  of  the  grid  points.  For  the  second  half  of  the 
k1"*1  iteration,  the  updating  of  field  values  is  reversed. 
Expressed  mathematically,  the  algorithm  is: 

f^ij1  =  fii j  +  '  nij)J  q  =  (5. 1.3. 1.1) 

[fl^j  =  ]  q  =  Z  (5. 1.3. 1.2) 

x  “  in 


where  l  *  -1,1  so  that  when: 
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Q  ~  " 1  j  -  2,4,6,...,  n-1  (5. 1.3. 1.3) 

Q  *  1  _>  j  -  1,3,5,...,  n  (5. 1.3. 1.4) 

In  equation  5. 1.3. 1.1,  S  is  the  acceleration  parameter  based 
on  the  dominant  eigenvalue  and  W^j  is  a  matrix  of  relaxation 
parameters  to  be  discussed  later.  The  variable  is  just 
a  dummy  variable  representing  either  i/k  ^  of  G^. .  The  promo¬ 
tion  factor  S  is  related  to  the  dominant  eigenvalue  by: 

s  =  - a -  (5.1. 3.1. 5) 

1  -  Uxl 

where  a  is  a  scaling  factor  so  that  S  is  bounded  since  in 
extreme  cases,  |  |  is  very  close  to  unity.  In  practice  S  is 

not  applied  on  every  increment  of  the  loop  parameter  k  but 
only  every  20  increments. 

Woo  found  that  a  range  of  10-15  for  S  to  give  about  the 
same  rate  of  convergence.  No  justification  was  given  for 
this  range.  It  appears  that  it  is  based  on  numerical  exper¬ 
imentation  since  no  eigenvalues  were  published  by  Woo.  In 
general  Woo  found  that  stability  is  not  sensitive  to  the 
convergence  promotion  factor  S  so  long  as  it  is  not  applied 
too  often.  When  the  convergence  promotion  factor  is  not 
being  applied,  the  factor  S  is  set  equal  to  one  in  Equation 
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The  use  of  relaxation  parameters  is  the  unique  features 
of  any  SOR  algorithm.  The  relaxation  factors  in  Equa¬ 
tion  5. 1.3. 1.1  plays  an  important  role  in  determining  the 
rate  of  convergence  of  the  solution.  An  optimum  value  usu¬ 
ally  exists,  however  for  most  cases,  there  exists  no  theoreti¬ 
cal  basis  for  calculating  this  optimum  value,  except  for  the 
particular  case  of  a  square  grid  system.  Varga4*  treats 
these  special  cases. 

j  2  26 

Both  Woo^  and  Pitter  et  al .  ,  used  a  constant  relaxa¬ 

tion  factor  for  the  stream  function,  i.e.,  Equation  4.3.1. 

They  based  it  on  the  work  of  Russell44,  who  studied  the  solu¬ 
tion  of  the  Navier- Stokes  equation  for  flow  over  flat  plates. 
Russell  suggested  that  the  relaxation  factor  for  the  stream 
function  equation  can  be  selected  as  uniform  over  the  entire 
field.  His  choice  of  relaxation  factor  is  based  on  the  num¬ 
ber  of  steps  of  each  independent  variable  and  for  the  case 
considered  here,  is  given  by: 


W.  = 


1  +  7T2(m'Z  +  n'2)1/2 


(5. 1.3. 1.6) 


where  m  and  n  have  been  previously  defined.  Pitter  tried  to 
use  Equation  5. 1.3. 1.6  for  the  case  of  flow  around  an  oblate 
spheroid  but  found  by  trial  and  error  that  using: 

W*  fT  +  TT(m'2  +  n'Z)1/Z 


(5. 1.3. 1.7) 
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gave  better  results.  Equation  5. 1.3. 1.7  is  used  in  the  cur¬ 
rent  code  NGCEFF  to  start  the  solution.  However,  after 
approximately  20  increments  of  the  loop  parameter  k,  is 
recomputed  using  the  following  relationship: 


W, 


opt 


i .  a  =)172 


(5. 1.3. 1.8) 


where  E  is  determined  as  follows:  Define  the  norm  of  vari¬ 
able  ip  as: 


i  k+1 


n  m 
z  Z 
j=l  i-1 


*k+1 


.k 

ip .  . 


Now  define  5  as: 


,  k  +  1 

s '  ™  A 


(5. 1.3. 1.9) 


(5.1.3.1.10) 


provided  ||  •  |j^  +  1  <  ||  •  ||*  and  ||  •  ||*  f1  0,  which 
indicates  convergence.  This  was  found  to  accelerate  con¬ 
vergence  of  the  ip  field  even  better  than  Equation  5. 1.3. 1.7. 

The  relaxation  of  the  vorticity  field  is  much  more 
difficult.  Equation  4.3.2  is  programmed  using  Equations 
5. 1.3. 9  -  5.1.3.11,  and  the  DEM  iteration  procedure  explained 
earlier.  However,  instead  of  having  a  constant  relaxation 
factor  for  the  entire  vorticity  field,  a  matrix  of  relaxation 
factors,  W-j ,  is  calculated  using  the  relationships  developed 
by  Woo^.  According  to  Woo,  Equation  4.3.2  could  be  reasonably 
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relaxed  to  a  numerical  solution  by  the  relationship: 


2  + 


/T"  Re  cosh  £ 
sinn  cosh£ 


(  cf*)2  * 


(5.1.3.1.11) 


i  =  2,3,...,  m-1 
j  —  2,3,...,  n- 1 


where  Re  is  the  Reynolds  number,  n  is  the  angular  coordinate 
and  £  is  the  radial  coordinate.  These  relaxation  factors 
were  used  but  because  of  severe  stability  problems  associated 
with  the  vorticity  transport  equation,  it  was  necessary  to 
stabilize  the  vorticity  field  while  the  stream  function  field 
was  allowed  to  relax.  Once  the  stream  function  field  con¬ 
verged,  the  vorticity  field  was  allowed  to  change  slowly.  Of 
course  this  caused  the  stream  function  field  to  change  which 
required,  because  of  stability  problems,  to  slow  down  the 
rate  of  convergence  of  the  vorticity  field.  Thus  the  matrix 
of  relaxation  factors  were  initially  multiplied  by  10  ** 
and  as  iteration  continued  the  multiplication  factor  was 
raised  until  it  reached  one. 

This  procedure  did  not  always  work  and  in  fact,  too 
often  the  iteration  procedure  would  diverge  rapidly  regard¬ 
less  of  the  value  of  relaxation  factors.  If  sensed  in  time 
the  divergence  could  be  stopped  by  multiplying  by  zero  and 
allowing  the  stream  function  field  to  converge  again,  and 


then  allow  small  changes  in  the  vorticity  field.  This  pro¬ 
cedure  rarely  worked  but  was  tried  before  stopping  the  com¬ 
puter  code. 

5. 1.3. 1.1  Properties  of  Iterative  Matrices  Being  Considered 

The  system  of  equations  generated  by  Equations  5. 1.3.1 
and  5. 1.3. 9  -  5.1.3.11  can  be  written  in  the  general  form: 

A  x  =  b  (5. 1.3. 1.1.1) 

or  converted  into  an  equivalent  form: 

x  =  T  x  +  c  (5. 1.3. 1.1. 2) 

where  A  is  the  matrix  of  coefficients,  b  is  a  column  vector 
of  boundary  conditions,  x  is  a  column  vector  of  the  unknowns, 
and  matrix  T  and  column  vector  c  are  quantities  which  depend 
on  which  matrix  iterative  procedure  is  used,  i.e.,  Jacobi, 
Gauss-Seidel ,  SOR,  etc. 

The  solution  to  Equation  5. 1.3. 1.1. 2  exists  and  is 
unique  provided  A  is  nonsingular.  For  small  systems  of  equa¬ 
tions,  its  solution  can  be  found  easily  by  computing  A'*  how¬ 
ever  for  large  systems  of  equations,  it  is  neither  practical 
nor  computationally  efficient  to  do  so. 

For  the  equations  considered  earlier  in  this  chapter, 
the  number  of  unknowns  and  hence  the  number  of  equations,  Nr, 


1 


Nr  =  (m  -  2) (n  -  2) 


which  makes  A  and  T  Nr  x  Nr  matrices. 


(5. 1.3. 1.1. 3) 


Matrix  A  in  Equation  5. 1.3. 1.1.1  can  be  expressed  as  a 


matrix  sum: 


A  =  D  -  L  -  U 


(5. 1.3. 1.1. 4) 


where  D,  L  and  U  are  the  diagonal,  lower  and  upper  triangular 
matrices.  In  matrix  notation,  the  point  Jacobi  interactive 


method  is  expressed  as 


Xk+1  =  B  xk  +  D"1  b 


where  B  is  the  point  Jacobi  matrix: 


B  =  D~ 1  (L  +  U) 


(5. 1.3. 1.1. 4) 


(5. 1.3. 1.1. 5) 


and  D"  is  the  inverse  of  D.  Even  though  the  point  Jacobi 
scheme  is  not  used  in  the  code  NGCEFF,  it  plays  an  important 
role  in  the  analysis  of  the  other  iterative  techniques. 

The  SOR  and  DEM  iterative  methods  in  matrix  notation 


xk  +  1  =  L  xk  +  (I  -  W  D'1  L)  D~  ^  W  b 


(5. 1.3. 1.1.6) 


where  L  is  the  point  successive  relaxation  matrix  given  by: 
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L  MI  ■  y'1  L)  ~ 1  {(I  -  W)I  +  W  D'1  U}  (5. 1.3. 1.1. 7) 

and  W  is  either  a  diagonal  matrix  corresponding  to  the  relaxa¬ 
tion  factors  multiplied  by  the  DEM  factor  or  a  scalar  value 
equal  to  the  relaxation  factor  multiplied  by  the  DEM  factor. 

Convergence  depends  on  whether  or  not  the  spectral  radius, 
p(A),  is  less  than  one.  To  each  of  the  iterative  methods 
described,  an  associated  error  vector  e  ,  defined  by: 

£k  -  xk  -  x  (5. 1.3. 1.1. 8) 

k  £  0 

where  x  is  the  solution  vector,  is  formulated.  Equation 
5. 1.3. 1.1. 8  can  be  rewritten  as: 

ek  -  Tk  c°  (5. 1.3. 1.1. 9) 

k  >  0 

where  T  is  the  matrix  B  or  L.  These  iterative  matrices  tend 
to  converge  or,  stated  another  way,  given  any: 

e°  >  0 


then 


ek  -*>  0,  if  p (M)  <  1.  (5.1.3.1.1.10) 

To  calculate  the  spectral  radius  of  the  point  successive 
relaxation  matrix  is  difficult  and  can  be  avoided  by  applying 
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the  following  theorem. 

Theorem:  Let  the  Jacobi  matrix  be  a  non-negative  nxn 

matrix  with  zero  diagonal  entries,  and  let  L  be  the  Gauss- 
Seidel  matrix,  the  special  case  W  =  (W)^j  =  1.  Then,  one 
and  only  one  of  the  following  mutually  exclusive  relations  is 
valid: 

1.  p (B)  =  p  (L)  =  0 

2.  0  <  p (L)  <  p (B)  <  1 

3.  1  =  p (B)  =  p(L) 

4.  1  <  p(B)  <  p (L) .  (5.1.3.1.1.11) 

The  above  theorem,  unfortunately,  does  not  apply  when 
matrix  B  is  not  a  non-negative  matrix.  Examples  can  be  for¬ 
mulated  where  one  iterative  method  is  convergent  while  the 
other  is  divergent,  when  B  is  not  a  non-negative  matrix. 
However,  if  a  matrix  is  not  non-negative  it  can  still  be  con¬ 
vergent  if  matrix  A  is  strictly  or  irreducibly  diagonally 
dominant . 

Besides  requiring,  p(A)  <  1,  a  measure  of  the  speed  or 
rate  of  convergence  is  the  condition  number,  K(A) ,  which  is 
very  important  in  the  study  of  rounding  error,  also.  In 
other  words,  when  K(A)  is  large,  small  relative  errors  in  A 
can  lead  to  large  errors  in  the  solution. 


Thus,  two  numbers,  the  spectral  radius  and  the  condition 
numbers,  are  computed  by  the  Code  NGCEFF  and  these  will  be 
reported  in  the  next  chapter. 

5.2  The  NGCEFF  Code 

The  numerical  determination  of  the  gravitational  colli¬ 
sion  shape  factor  0  is  accomplished  by  the  NGCEFF  code,  which 
calculates  both  the  nonspherical  and  its  volume  equivalent 
spherical  gravitational  collision  efficiencies.  Figure  5.2.1 
is  a  flowchart  of  the  code.  The  code  first  calculates  the 
nonspherical  collision  efficiency,  e^,  based  on  the  Reynolds 
number  of  the  oblate  spheroid.  The  volume  equivalent  radius, 
based  on  the  semi-major  axis  length  of  the  oblate  spheroid, 
is  determined  and  a  Reynolds  number  is  guessed  for  a  spheri¬ 
cal  particle.  From  the  results  of  solving  the  Navier-Stokes 
equation,  a  comparison  is  made  between  the  calculated  size 
(radius),  based  on  the  guessed  Reynolds  number,  and  the 
actual  volume  equivalent  radius  for  the  sphere.  A  new 
Reynolds  number  is  guessed  if  the  size  and  actual  radius  dis¬ 
agree.  When  agreement  is  achieved,  the  spherical  gravita¬ 
tional  collision  efficiency.  Eg,  is  calculated  and  then  the 
shape  factor  0.  This  procedure  is  done  so  as  to  insure 
uniformity,  i.e.,  the  same  assumptions  and  algorithms  are 
used  to  determine  and  eg. 


The  NGCEFF  code  is  written  in  IBM  double  precision 
FORTRAN  language.  The  main  program  controls  the  logic  as 
depicted  in  Figure  5.2.1.  The  listing  of  the  code  is  given 
in  an  appendix.  Not  listed  are  subroutines  that  are  part 
of  IMSL"^  and  UNPACK^  subprogram  libraries. 

Except  for  IMSL  and  UNPACK  subroutines,  this  section 
describes  the  subroutines  and  function  subprograms  of  the 
NGCEFF  code. 

5.2.1  Code  Structure 

Since  numerical  solutions  to  the  Navier-Stokes  equation 
requires  large  arrays,  the  code  uses  many  common  blocks  and 
equivalence  statements  so  that  arrays  can  be  used  as  working 
or  multi- variable  arrays.  This  is  accomplished  by  storing 
intermediate  results  that  are  needed  later,  on  direct  access 
storage  devices  such  as  disks.  Most  of  the  output  and  input 
from  disks  is  done  using  unformatted  I/O  statements,  which 
are  extremely  fast  data  transfer  statements  since  there  is 
a  one-to-one  correspondence  between  internal  storage  loca¬ 
tions  (bytes)  and  external  record  positions. 

Table  5. 2. 1.1  lists  the  common  blocks  and  gives  a 
description  of  the  type  of  variables  contained  in  it.  Since 
many  of  the  subroutines  are  common  to  both  spheroid  and 
spherical  particles,  a  labeled  common  block  SWITCH  is  placed 
in  those  subroutines  and  the  values  of  the  variables  are 
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changed  for  the  type  of  particle  being  considered.  Labeled 
common  blocks  OBLAT$  and  WORK  are  used  to  store  subscripted 
variables.  With  the  exception  of  SWITCH,  OBLAT$  and  WORK, 
the  values  of  the  variables  in  the  remaining  common  blocks 
do  not  change  once  determined  for  a  specific  case. 

5. 2. 1.1  Main  Program 

The  MAIN  program  reads  in  the  values  of  variables  needed 
to  run  the  code  and  controls  its  operation.  Before  termina¬ 
ting  the  code,  MAIN  checks  to  determine  if  another  case  needs 
to  be  evaluated.  Table  5. 2. 1.1.1  lists  the  input  variables. 
Variable  GOTO  allows  the  user  to  consider  very  similar  cases 
without  having  to  redefine  all  the  input  variables.  The  MAIN 
program  calls  subroutine  INITAL,  OBLATE,  SIZE,  control  parame¬ 
ters.  The  first,  I,  identifies  the  type  of  large  particle  be¬ 
ing  considered  and  how  much  information  is  known  about  that 
particle.  The  second  parameter,  J,  is  used  to  indicate  that 
only  the  variables  associated  with  the  similar  particle  have 
been  changed.  This  allows  certain  subroutines  to  be  skipped. 
The  final  parameter,  JSTART,  is  used  to  generate  the  initial 
guess  for  the  stream  function  and  vorticity  fields,  and  then 
is  used  to  check  that  the  radius,  drag  coefficient,  and 
Reynolds  number  all  agree  for  the  volume  equivalent  particle. 
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Table  5. 2. 1.1.1 


Input  Variables  to  NGCEFF  Program 


000  to 

00020 

c 

00030 

c 

00040 

c  input: 

PARAMETER 

CARD ( COLUMNS > 

FORMAT 

USAGE 

oooso 

c 

MU 

1 

<1-10 

) 

010.0 

AIR  VISCOSITY 

00060 

c 

RHG 

1 

<11-20 

) 

010.0 

AIR  DENSITY 

00070 

c 

LAMBDA 

1 

<21.-30 

> 

010.0 

MEAN  FREE  PATH  OF  AIR 

00080 

c 

EPS 

1 

<31-40 

> 

<510.0 

ERROR  CRITERIA  FOR  Al.l. 

00090 

C****» 

00100 

c 

REYN1 

2 

<1-10 

) 

010.0 

REYNOLDS  NUMBER  FOR  OBLATE 

oouo 

c 

AR1 

2 

<1-10 

) 

GtO.O 

AXIS  RATIO  OF  OBLATE 

00120 

c 

DENS1 

2 

<11-20 

) 

G10.0 

BULK  DENSITY  OF  OBLATE. 

00130 

c 

AI.PHA1 

2 

<21-30 

) 

G10.0 

DENSITY  CORRECTION  FACTOR 

00140 

c 

DISK 

2 

<31-40 

) 

1510.0 

SELECT  OPTION  NEEDED... 

00  ISO 

c 

0 . 0-OEHERA TE  PSI  AND  0  FIELDS 

00160 

c 

1.0-READ  FIELDS  FROM  DISK. 

00170 

c 

UNIT  NUMBER  IS  FT09F001 

00160 

cat*** 

00190 

c 

DZETA 

3 

<1-10 

) 

G10.0 

RADIAL  FlTEP  SIZE.  SPHEROIDAL 

00200 

c 

MP1 

3 

<11-15 

> 

IS 

NUMBER  OF  ANGULAR  STEPS 

00210 

c 

N 

3 

<16-20 

> 

15 

NUMBER  OF  RADIAL  STEPS 

00220 

c 

MATRX 

3 

<21-25 

> 

15 

SELECT  OPTION  NEEDED... 

00230 

c 

O-NO  MATRIX  ANALYSIS 

00240 

c 

1-HAIRIX  analysis: 

00230 

c 

SPECTRAL  RADI  LIS  .EIGENVALUES 

00240 

c 

AND  KECIPOCAL  CONDITION 

00270 

c 

NUMBER  CALCULATED 

00280 

c 

IOUT 

3 

<26-30 

) 

15 

SELECT  OPTION  NEEDED... 

00290 

c 

1-OUTPUT  RESULTS 

00300 

c 

FT08F001  MUST  BE  DEFINED 

00310 

c 

FT09F001  MUST  »F:  DEFINED 

00320 

c 

FT  1  OF 001  MUSI  BE  DEFINED 

00330 

f; 

ABOVE  DATA  SETS  USE 

00340 

c 

UNFQRHATF.D  WRITE  «  READ 

00350 

c 

STATEMENTS 

00360 

c 

O-NO  OUTPUT  NEEDED 

00370 

etc*** 

00380 

c 

R2V 

4 

<1-10 

> 

r,io.o 

VOLUME  EDUI VALENT  RADIUS  OF 

00390 

c 

SMALLER  PARTICLE. 

00400 

c 

dens;,1 

4 

<11-20 

> 

010.0 

BULK  DENSITY  OF  SMALLER 

00410 

c 

AI.PHA2 

4 

<21-30 

) 

G10.0 

DENSITY  CORRECTION  FACTOR  OF 

00420 

c****» 

00430 

c 

GOTO 

5 

<1-10 

) 

G10.0 

SELECT  OPTION... 

00440 

c 

1.0-NEW  CARDS  1-4  NEEDED. 

00450 

r. 

2.0-NEU  CARDS  2-4  HEEDED. 

00460 

c 

3.0-MEU  CARDS  3*4  NEEDED 

00470 

c 

4.0-0NLY  CARD  4  CHANGED 

00480 

c 

99.0-STOP  CODE. 

00490 

c 

00500 

c 

00510 

c*********************************************************************** 

Ill 


5. 2. 1.2  Subroutine  INITAL 

*  The  subroutine  INITAL  initializes  the  parameters  neces¬ 

sary  to  carry  out  calculations  by  the  subroutine  COLL.  Not 
all  parameters  are  determined  the  first  time  the  subroutine 
is  called,  since  not  all  the  information  is  known  and  there¬ 
fore  subroutine  INITAL  is  called  several  times  as  more  infor¬ 
mation  is  determined.  The  following  is  a  list  of  the  para¬ 
meters  and  their  definitions: 


GV 

-  Gravitational  constant 

PI 

-  The  mathematic  constant  it 

CLL 

-  Collision  parameter 

ERR 

-  Error  criterion  for  bisection  method 

in 

COLL 

CDTOTL 

-  Drag  coefficient  OBLATE 

NEQTN 

-  Number  of  equations,  IMSL  subroutine 

DVOGER 

MAXDER 

-  Maximum  order,  IMSL  subroutine  DVOGER 

HD 

-  Initial  time  step  size,  IMSL  subroutine 

DVOGER 

HDMIN 

-  Minimum  time  step  size,  IMSL  subroutine 

DVOGER 

HDMAX 

-  Maximum  time  step  size,  IMSL  subroutine 

DVOGER 

MTH 

-  Integration  method,  IMSL  subroutine 

DVOGER 

GAMMA 

-  Ratio  of  smaller  particle  density  to 

larger 

particle  density 

BRHO 

-  Ratio  of  particle  buoyancy  effects 

IAR1 

-  Ratio  of  inertia  effects  for  oblate 

spheroid 

particle 

I 
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IR1V 

VINF2 

REYN2 

CDAR1 

SEMI 

R1V 

UN  IF 

VINF1 

R 

STK1 
LYK1 
CDV1 
REYN VI 
VINFV1 

STKV1 

LYKV1 


Ratio  of  inertia  effects  for  volume  equivalent 
particle 

Stokes  terminal  velocity  for  smaller  particle 
Reynolds  number  of  smaller  particle 
Drag  coefficient  for  smaller  particle 
Semi-major  axis  length  for  oblate  spheroid 
Radius  of  volume  equivalent  particle 
Stokes  terminal  velocity  of  larger  particle 
Actual  terminal  velocity  of  oblate  spheroid 
Ratio  of  volume  equivalent  radii 
Stokes  number  based  on  oblate  spheroid 
Interaction  number  based  on  oblate  spheroid 
Drag  coefficient  of  volume  equivalent  particle 
Reynolds  number  of  volume  equivalent  particle 
Actual  terminal  velocity  of  volume  equivalent 
particle 

Stokes  number  based  on  volume  equivalent 
particle 

Interaction  number  based  on  volume  equivalent 
particle . 


5. 2. 1.3  Subroutine  OBLATE 


As  the  name  implies,  this  subprogram  solves  the  Navier 
Stokes  equation  for  flow  around  an  oblate  spheroid.  Since 
an  axis  ratio  of  0.999  has  been  shown  to  approximate  a 
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sphere  ,  this  subroutine  also  gives  the  solution  for  spheres. 
The  equations  used  to  find  the  stream  function  and  modified 
vorticity  fields  are  given  in  Section  5.1.3. 

To  start  the  iteration  scheme,  which  has  been  described 
in  the  previous  chapter,  an  initial  guess  is  made  for  the 
stream  function  (PSI)  and  the  modified  voriticity  (G)  fields. 
Two  options  are  available.  If  no  previous  solutions  are 
available,  then  DISK  should  equal  0  so  that  subroutine  START 
is  called  to  generate  an  initial  guess  by  using  Oberbeck's 
solution.  If  DISK  does  not  equal  zero,  the  subroutine  OBLATE 
assumes  that  starting  (initial  guess)  values  are  on  a  random 
access  device  with  a  file  name  FT09F001.  When  convergence 
is  achieved,  control  is  returned  to  the  MAIN  program  with 
J START  equal  to  2.  Key  parameters  returned  are  PSI,  G,  and 
the  total  drag  coefficients,  CDTOTL. 

The  pressure  drag  coefficient  Cp  is  determined  from  the 
relation : 

TT 

Cp  =  2  /  k  sinri  cosq  dn  (5. 2. 1.3.1) 

o 

and  the  friction  drag  coefficient  is  determined  by 

Cf  -  ^  tanh£0  GQ  sinq  dp  (5. 2. 1.3. 2) 

where  the  total  drag  coefficient  CDTOTL  is  just  the  sum  of 
the  pressure  drag  and  friction  drag  coefficients.  The 
terms  in  Equations  5.2.1. 3.1  and  5.2.1. 3.2  are 


k  ~  Nondimens ional  pressure  unit 


i  +  &  f  b  Boj  I 

Ee  j.  3n  I  n=0  ^ 


■  3  to 


•  [-5-^  +  (jjtanh^Jdn 

3C  C  -C. 


[5. 2. 1.3. 3) 


Gq  ~  Modified  vorticity  at  surface  oblate  spheroid, 
i  .  e .  ,  £  =  £ 


10  ~  Vorticity 


Re  ~  Reynolds  number 

~  Value  of  5  at  the  outer  envelope. 

Equation  5. 2. 1.3. 3  was  evaluated  by  the  trapezoidal  rule, 
since  higher  order  integration  methods  were  unstable. 

The  iteration  scheme  in  subroutine  OBLATE  for  finding 
the  stream  and  vorticity  fields  was  not  as  stable  as  was 
hoped.  It  is  very  sensitive  to  the  relaxation  parameters 
and  the  initial  guess.  Thus  if  the  Reynolds  number  is  not 
approximately  equal  to  0,  the  use  of  subroutine  START  is 
not  recommended. 
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5. 2. 1.4  Subroutine  SIZE 

This  subroutine  is  used  only  to  determine  the  correct 
Reynolds  number  for  the  volume  equivalent  particle.  Once 
the  semi-major  axis  of  the  spheroid  is  determined,  the 
volume  equivalent  radius  is  known  but  not  the  Reynolds  num¬ 
ber.  As  an  initial  guess,  Stokes  law  can  be  used  and  then 
subroutine  OBLATE  is  called.  The  results  of  OBLATE  allow 
the  radius  of  a  sphere  corresponding  to  the  guessed  Reynolds 
number  to  be  calculated.  Comparison  of  this  radius  with  the 
volume  equivalent  radius  allows  the  program  to  adjust  the 
Reynolds  number.  Subroutine  OBLATE  is  then  called  again. 
After  three  guesses,  subroutine  SIZE  uses  ISML  subroutines 
ICSICU  and  ICSEVU  to  interpolate  the  next  Reynolds  number. 
This  procedure  is  continued  until  the  absolute  difference 
between  the  calculated  radius  and  the  volume  equivalent 
radius  is  less  than  EPS.  When  this  occurs,  control  JSTART 
is  set  equal  to  four  and  the  two  arrays  PSI  and  G  are  placed 
on  a  random  access  device  using  file  name  FT09F001.  This  is 
done  so  that,  if  the  user  wants  to  consider  similar  cases, 
the  initial  guess  can  be  closer  to  the  actual  answer  than 
the  analytically  generated  guess. 

5. 2.1. 5  Subroutine  VELCTY 

Given  the  discrete  values  of  the  stream  function  in 
spheroidal  coordinates,  this  subroutine  will  calculate  the 
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velocity  fields  in  cylindrical  coordinates  by  using  Stokes 
relationships  and  the  transformation  expressions  listed  in 
Chapter  4.  Results  are  stored  on  a  random  access  device 
using  file  name  FT10F001. 

5. 2. 1.6  Subroutine  COLL 

The  COLL  subroutine  utilizes  the  results  from  INITAL 
and  VELCTY  and  integrates  the  six  nondimens ional  equations 
of  aerosol  motion  (Equations  4.4.1.34  -  4.4.1.39)  to  deter¬ 
mine  eg  and  e^.  The  velocity  fields,  which  are  necessary 
for  estimating  drag  forces,  are  read  in  from  a  random  access 
device  using  file  name  FT10F001.  The  older  IMSL  routine 
DVOGER  is  the  integrator  because  it  offers  some  advantages 
over  the  newer  integrator  routine  DGEAR. 

The  initial  conditions,  Equations  5. 1.1. 7  -  5.1.1.1.11, 
are  computed  and  control  parameters  are  initialized.  DVOGER 
is  repeatedly  called  to  determine  the  particles'  trajectory 
for  an  initial  horizontal  separation.  Integration  is  con¬ 
tinued  until  a  hit  or  miss  occurs  between  the  larger  parti¬ 
cle  and  the  smaller  particle.  COLL  determines  the  angle  of 
separation  and  distance  of  separation  (Equations  4.2.2.10  - 
4.2.2.13)  for  each  initial  horizontal  separation.  Whenever 
a  grazing  collision  occurs,  i.e.,  distance  of  separation  is 
less  than  10-6,  the  horizontal  separation  distance  is  saved. 
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As  indicated  in  Chapter  III,  after  the  maximum  critical 
grazing  trajectory  is  found,  it  is  still  necessary  to  deter¬ 
mine  the  inner  critical  trajectory  offset.  Once  this  is 
accomplished,  the  gravitational  collision  efficiency  can  be 
computed. 

The  search  strategy  is  to  start  with: 


p)Q  =  (a*  +  *1>v)/r*>y  (OBLATE) 


(5. 2. 1.6.1) 


p)o  =  (r^  v  +  r2  v)/r*i  v  (Volume  bquivalent  (5. 2. 1.6. 2) 

Sphere) 

where  p)Q  is  the  nondimensional  initial  horizontal  offset, 
a  is  the  semi-major  axis  length,  r,  is  the  large  sphere 

1  ,  V 

volume  equivalent  radius,  and  r7  is  the  volume  equivalent 

L  y\J 

radius  of  the  smaller  sphere.  To  continue  searching,  the 

Y 

new  horizontal  offset,  p)Q,  is: 


p)o  _  p)o  (1  '  TTT) 


(5. 2. 1.6. 3) 


k  =  1,2 .  10 


until  a  hit  occurs.  Once  this  happens  the  bisection  tech¬ 
nique  is  used  to  determine  the  bisection  horizontal  offset, 
p) b ,  which  is : 


P)b  ■  [P)^1  *  P)*]/2 


(5. 2. 1.6. 4) 
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The  search  strategy  for  the  minimum  horizontal  offset  is  the 
same  as  for  the  maximum  horizontal  offset. 

The  advantages  of  the  DVOGER  routine  is  that  it  will 
adjust  the  time  step  size  internally  to  meet  error  criteria 
and  to  achieve  finer  numerical  detail  of  the  particle  tra¬ 
jectory  calculations.  In  addition,  DVOGER  has  the  capa¬ 
bility  to  repeat  the  last  step  of  integration  if  for  some 
reason  too  large  of  a  time  step  was  taken  which  caused  the 
larger  particle  to  pass  through  the  smaller  particle  (physi¬ 
cally  not  realistic).  When  these  problems  occur,  new  maxi¬ 
mum  and  minimum  step  sizes  can  be  given  to  DVOGER  in  order 
to  complete  the  integration.  These  features  are  not  avail¬ 
able  in  DGEAR. 

5. 2. 1.7  Subroutine  OUTPUT 

The  OUTPUT  subroutine  summarizes  all  the  input  data  and 
the  key  intermediate  and  final  results.  It  is  called  three 
times  by  the  MAIN  program.  The  routine  first  prints  the 
values  of  the  input  parameters  associated  with  control  of 
the  code  and  those  quantities  for  containment  conditions. 
Later  OUTPUT  is  called  to  give  the  nonspherical  gravitational 
collision  efficiency  eN  and  particle  quantities  used  to  cal¬ 
culate  the  efficiency.  The  last  time  OUTPUT  is  called 
results  in  giving  spherical  gravitational  collision  effic¬ 
iency  £g  and  the  collision  shape  factor  6. 
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5. 2. 1.8  Subroutine  IDERIV 


Subroutine  IDERIV  is  called  twice  by  subroutine  VELCTY 
to  calculate  the  derivatives  of  the  stream  function  field 
(PSI)  .  IDERIV  uses  the  cubic  spline  coefficients  determined 
by  IMSL  routine  ICSICU.  The  two  derivatives  determined  by 
IDERIV  are  those  in  Equation  4.3.7  and  4.3.8.  Basically 
IDERIV  is  the  IMSL  routine  ICSEVU  except  the  derivative  of 
a  cubic  spline  is  calculated  instead  of  the  IMSL  library 
notebook,  the  value  of  the  derivative  of  the  spline  approxi¬ 
mately  at  the  points  IK  is: 

Sr  =  (3C.  ,D  +  2C -  7)D  +  C.  ,  (5. 2. 1.8.1) 

1  3  c3  J  > L  J  *  •*- 


where : 


S7  -  Derivatives  at  the  points  IK  ,  i  =  l,2,...  m 


C .  ,  ~  Spline  coefficients,  k=l,2,3;  j  =  l,2,...,  nx-1 
J  * K 

nx  ~  Number  of  function  values  of  data  point 
m  ~  Number  of  interpolating  points 
and  D  =  IK  -  Xj  where  the  interval  is  determined  by: 

Xj ;  ui  ^  Vi 

where  Xj  are  the  abscissae  data  points.  If  the  derivatives 
are  needed  only  at  the  original  data  points,  then  D  =  0  and 
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5.2. 1.9  Subroutine  DEFINE 

This  subroutine  loads  the  proper  values  into  the  common 
block  called  SWITCH.  This  allows  subroutines  like  OBLATE, 
VELCTY  and  COLL  to  be  used  by  either  the  oblate  spheroidal 
particle  or  its  volume  equivalent  spherical  particle.  Para¬ 
meters  associated  with  the  oblate  spheroid  are  in  common 
block  AERSL1  and  those  associated  with  the  volume  equivalent 
sphere  are  in  common  block  VER1.  Subroutine  DEFINE  uses  the 
MAIN  control  parameter  I  to  correctly  load  SWITCH. 

5.2.1.10  Subroutine  DFUN 

The  subroutine  DVOGER  requires  a  user  supplied  sub¬ 
routine  which  contains  the  differential  equations  to  be 
integrated.  Thus  DFUN  contains  the  six  nondimensional  equa¬ 
tions  of  aerosol  particle  motion,  Equations  4.4.1.34  - 
4.4.1.39. 

The  drag  force  terms  used  in  the  equation  are  calculated 
in  separate  subroutines.  Since  the  explicit  drag  force  terms 
are  not  given  in  DFUN,  but  are  calculated  elsewhere,  this 
subroutine  is  very  general  and  allows  modification  of  the 
drag  force  terms  in  the  separate  subroutines  without  affect¬ 
ing  DFUN. 

The  DVOGER  subroutine  allows  the  user  to  provide  a  sepa¬ 
rate  subroutine  to  calculate  the  Jacobian  of  the  equations  of 
motion.  This  occurs  when  MTH  is  equal  to  1.  Since  analytical 
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expressions  are  not  available  for  the  drag  force  terms,  the 
Jacobian  has  to  be  calculated  numerically.  This  is  automat¬ 
ically  done  by  the  subroutine  DVOGER  if  MTH  is  equal  to  2. 
The  array  associated  with  the  Jacobian  is  PW(6,6)  and  can  be 
used  as  a  working  array. 

Drag  forces  are  called  by  subroutine  FORSUP,  which 
utilizes  function  subprograms  WRHOL,  WZ1,  WRH02,  and  WZ2. 

5.2.1.11  Subroutine  FORSUP 

This  subroutine  is  called  by  subroutine  DFUN  to  calcu¬ 
late  the  drag  forces.  Drag  forces  are  determined  by  using 
Equations  1.4.2.1  and  4.4. 2. 2,  which  is  based  on  the  dif¬ 
ferences  between  the  nondimensional  particle  velocities  vpl» 
Vp2 ,  vzi>  VZ2»  an(*  t^ie  generated  nondimensional  velocity 

fields  W  ,,  W  ,  W  ,,  W  ,  which  are  determined  by  function 
pi ’  p2 ’  zl ’  zZ ' 

subprograms  WRH01,  WZ1,  WRH02,  WZ2. 

5.2.1.12  Function  Subprograms  WRH01,  WRH02,  WZ1,  WZ2 

These  subprograms  determine  the  velocity  around  the 
particle  at  the  point  in  space  where  the  other  particle  is. 
The  subprograms  WRH01  and  WZ1  use  the  results  from  VELCTY 
and  the  interpolating  Equations  5.1.2.2.16  and  5.1.2.2.17 
to  determine  and  Wpl,  respectively.  The  weighing  fac¬ 
tors  are  determined  by  the  formulas  in  the  same  section. 


122 


The  subprogram  WRH02  and  WZ2  use  Equations  4.4. 2. 3 
and  4.4. 2.4  to  calculate  Wj2  and  Wp2,  respectively. 

5.2.1.13  Function  Subprograms  STREAM  and  VORTCY 

These  subprograms  are  called  by  subroutine  START  to 
calculate  the  stream  function  and  vorticity  fields  around 
oblate  spheroids  and  spheres.  Oberbeck's  and  Stokes  solu¬ 
tions  are  imbedded  in  them. 


CHAPTER  VI 


RESULTS  AND  DISCUSSIONS 


6 . 0  Introduction 

This  chapter  utilizes  the  theory  developed  in  Chapter 
IV  and  the  numerical  methods  described  in  Chapter  V  to  cal¬ 
culate  values  of  the  gravitational  collision  efficiency  for 
spherical  and  nonspherical  particles  and  the  collision  shape 
factor  3. 

The  aerosol  sizes  and  densities  were  determined  by 
reviewing  the  theoretical  analysis  of  the  characteristics 
of  typical  LMFBR  particles  and  agglomerates.  No  attempt 
was  made  to  model  the  collisions  between  raindrops  and  ice 
particles  per  se ,  a  subject  of  great  interest  in  atmospheric 
science.  Fluid  parameters  were  selected  based  on  predicted 
HCDA  post-accident  conditions  inside  the  LMFBR  containment. 

Some  important  results  presented  in  this  chapter  are 
(1)  the  instability  of  the  numerical  method  used  to  solve 
the  vorticity  transport  equation  and  how  to  deal  with  this 
problem,  (2)  the  effect  caused  by  incorrectly  determining 
the  separation  distance  between  the  spherical  and  nonspher¬ 
ical  particles  before  collision  occurs,  and  (3)  the  results 
of  the  program  NGCEFF,  i.e.,  eg,  e^,  3  and  other  key  quanti¬ 
ties  are  reported. 
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6 . 1  Some  Important  Numerical  Methods  Results 

6.1.1  Solution  to  the  Navier-Stokes  Equation 

Recalling  Equations  4.3.1  and  4.3.2  and  the  numerical 
techniques  discussed  in  Section  5.1.3,  solutions  to  the 
Navier-Stokes  equations  were  sought,  starting  with  Re  =  0 
and  then  increasing  it  to  Re  =  5. 

Convergence  was  fairly  rapid  when  Oberbeck's  solution 
was  used  as  the  trial  solution  and  the  Re  =  0 .  This  occurs 
because  the  vorticity  transport  equation  is  reduced  to: 

E2G  =  0  (6. 1.1.1) 

and  therefore  the  Jacobian  operator  does  not  need  to  be 
computed.  The  stream  function  field  converged  fairly 
rapidly  (less  than  200  iterations) ,  but  the  vorticity  field 
was  more  difficult  to  relax.  Relaxation  factors  in  Equation 
5.1.3.1.11  were  multiplied  by  0.001  initially  to  prevent 
divergence.  Once  a  stable,  non-oscillating  field  was  achieved, 
the  relaxation  scaling  factor  could  be  increased  until  it  was 
greater  than  0.7. 

Although  a  converged  solution  was  achieved  for  Ober¬ 
beck's  solution,  the  use  of  this  solution  for  Reynolds  num¬ 
bers  greater  than  zero  proved  to  be  impossible.  Convergence 
was  more  dependent  on  trying  different  combinations  of  scal¬ 
ing  factors  and  DEM4^  factors  than  on  the  internal  logic  of 


the  numerical  method.  It  was  decided  that  further  investi¬ 


gations  were  necessary  than  reported  by  other  investiga¬ 
tors26,42. 

Before  proceeding  it  is  worthwhile  to  discuss  the  moti¬ 
vation  behind  this  investigation.  In  calculating  the  gravi¬ 
tational  collision  efficiency  the  subroutine  OBLATE  is  used 
to  calculate  the  drag  coefficients  for  both  the  nonspherical 
particle  and  the  volume  equivalent  sphere.  Furthermore,  the 
determination  of  the  size  of  these  two  particles  required  an 
inner  and  outer  iteration  scheme.  The  outer  iteration  scheme 
was  between  the  subroutines  SIZE  and  OBLATE  to  adjust  the 
Reynolds  number,  which  is  not  known  a  priori ,  to  the  drag 
coefficient  from  the  solution  of  the  Navier-Stokes  equation. 
When  these  two  quantities  are  known,  the  volume  equivalent 
radius  of  the  oblate  spheroid  can  be  compared  with  the  one 
determined  from  Equation  5. 1.1. 1.5  and 

2  2  1/3 

rl,v  =  {3p  CD  Rev,l  AR(P1  '  p)  pg/32}  ’  (6. 1.1. 2) 


The  Reynolds  number  must  be  changed  until  the  three  quanti¬ 
ties  (Rey  CD,  r^  v)  give  consistent  solutions.  The  inner 
iteration  procedure  is  the  Navier-Stokes  numerical  scheme, 
which  typically  required  between  1500-2000  iterations  under 
the  best  conditions  before  convergence  was  achieved.  Thus 
a  numerical  method  had  to  be  choosen  which  could  be  pre¬ 
dicted  to  give  convergence  without  user  intervention. 
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An  investigation  of  the  problem  proceeded  as  follows. 

A  review  of  the  stream  function  and  vorticity  field  lattice 
values  was  done  for  the  cases  where  the  stream  function  \J> 
was  held  constant  and  vorticity  allowed  to  vary  until  con¬ 
vergence  or  divergence  occurred.  This  procedure  was  repeated 
for  different  values  of  Reynolds  numbers.  A  similar  techni¬ 
que  was  used  for  the  stream  function  field,  holding  the 
vorticity  field  constant.  Actual  convergence  was  possible 
when  the  vorticity  field  was  held  constant;  however,  allow¬ 
ing  the  vorticity  field  to  change  after  the  stream  function 
field  had  converged  caused  divergence,  resulting  in  overflow 
problems  almost  immediately  if  the  size  of  the  Reynolds  num¬ 
ber  exceeded  0.001. 

At  this  point  the  spectral  radius  and  condition  number 
of  the  vorticity  transport  matrix  was  determined  for  a  20 
x  20  matrix  (n  =  6,  m  =  7,  see  Equation  5 . 1 . 3 . 1 . 1 . 3)  . 

Results  are  given  below: 


Re 

p(oj) 

KO) 

pO) 

KO) 

0.5 

1431 

4  x  10'7 

4.199 

1.2  x 

10 

0.0001 

4.48 

1.2  x  10'1 

4.199 

1.2  x 

10 

where  the  AR  =  0.999  and  AC  =  1.6.  The  problem  of  conver¬ 
gence  is  immediately  obvious  for  Reynolds  number  greater 
than  0.001. 
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Since  these  matrices  do  not  have  spectral  radii  less 
than  one  (see  Chapter  V)  any  oscillates  in  the  numerical 
technique  applied  will  cause  problems. 

A  close  review  of  the  voriticity  field  indicated  that 
the  oscillations  started  near  the  boundary  of  the  oblate 
spheroid.  The  boundary  condition  at  the  surface  of  the 
spheroid  is. 


a)  =  G/sinn,  5  =  iQ  (6. 1.1. 3) 

which  is  evaluated  by  means  of  Equation  4.3.1,  and  when 
written  in  finite  difference  form  becomes 


G 


i ,  1 


cosh  E 

o 

- 2 - 2 

sinh  E0  +  cosp 


(6. 1.1. 4) 


2  2 

It  is  the  representation  of  9  i p/dE,  on  the  surface  of  the 

spheroid  which  causes  oscillations  of  the  vorticity  field, 

which  if  not  controlled  by  methods  used  by  Woo4^  and  Pitter 
2  6 

et  al .  ,  will  lead  to  divergence.  Pitter  et  al .  used  the 

following  one-sided  finite-difference  equation, 


8  *i,l 


8  <1^ 


2  "^i  3 

0 


2(H)' 


(h) 


which  when  replaced  by  the  expression; 


(6. 1.1. 5) 
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did  not  cause  oscillates  near  the  boundary.  Higher  order 
forward  representations  were  derived  by  replacing  succes¬ 
sively  more  terms  in  the  Taylor  series  expansions  and  making 
use  of  the  fact  that  ip-  .  =  0  on  the  boundary.  These 

1  9  -L 

representations  caused  more  problems  by  introducing  higher 
order  harmonics  which  could  not  be  damped  out  by  the  scaling 
factors.  Thus,  Equation  6. 1.1. 6  is  the  expression  used  in 
NGCEFF . 


6.1.2  Solution  to  the  Dynamic  Equations  of  Motion  and 
Critical  Grazing  Path 

!  TO 

The  use  of  Gear  s  method  to  solve  the  six  equations 
of  particle  motion,  Equations  4.4.1.34  -  4.4.1.39,  was  sat¬ 
isfactory.  Remembering  that  the  equations  are  integrated 
for  the  purpose  of  determining  if  a  direct,  grazing,  or 
missed  collision  occurs,  it  is  necessary  to  start  the  inte¬ 
gration  routine  when  the  particles  are  at  a  distance  of 
at  least  50  radii  from  each  other.  When  this  was  done,  it 
took  approximately  200  steps  to  integrate  these  equations 
to  get  the  separation  angle  and  particle-particle  separation 
distance  (see  Figure  5. 1.1. 3).  Normally  fifteen  to  twenty 
initial  horizontal  offset  values  were  needed  before  the 
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grazing  path  was  determined  from  the  cubic  spline  curve 

fitting  technique  described  in  Section  5. 1.1. 1.2. 

One  key  requirement  in  solving  these  six  equations  of 

motion  is  calculation  of  the  forces  acting  on  the  particles. 

Interpolation  of  the  velocity  fields  must  be  done  correctly 

or  else  the  forces  will  be  incorrectly  applied  by  the  super- 

2  6 

position  principle.  Pitter  et  al .  incorrectly  programmed 
Equations  5. 1.2. 2.1  -  S. 1.2. 2. 2  letting  c,  the  characteristic 
length  in  the  spheroid  coordinate  system,  be  equal  to  cosh^Q, 
where  CQ  is  the  surface  of  the  oblate  spheroid.  This  problem 
was  first  determined  when  his  collision  trajectory  consis¬ 
tently  ended  up  inside  the  oblate  spheroid.  Time  did  not 
permit  modification  of  his  program  to  check  his  reported 
observation  of  an  annular  collision  domain,  but  it  is  sus¬ 
pected  that  this  programming  error  contributed  to  his  find¬ 
ings, 

6 . 2  Verification  of  NGCEFF  Routines 

It  was  important  to  verify  the  various  routines  in  the 
NGCEFF  code.  Each  major  subroutine  was  subjected  to  exten¬ 
sive  testing  before  merging  it  with  NGCEFF.  For  example, 
the  subroutines  VORTCY  and  STREAM  were  developed  to  generate 
the  vorticity  field  and  stream  function  values  for  Oberbeck's 
solution  (Re  =  0),  Verification  of  these  two  subroutines 
was  checked  by  comparing  the  generated  fields  with  those 
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from  solving  the  Navier-Stokes  equation  with  the  Reynolds 
number  equal  to  zero.  This  procedure  also  verified  the 
output  from  subroutine  OBLATE  since  the  initial  fields  used 
by  OBLATE  were  not  converged  fields. 

Subroutine  VELCTY  was  checked  by  taking  the  output 
from  OBLATE  and  STREAM  and  letting  the  subroutine  VELCTY 
operate  on  it  to  create  the  velocity  fields  in  oblate 
spheroidal  coordinate  system.  These  results  were  then  com¬ 
pared  to  the  analytical  solution  available  from  Oberbeck's 
solution  for  Reynolds  number  equal  to  zero.  This  procedure 
in  essence  also  verified  the  correctness  of  the  subroutines 
OBLATE  and  STREAM. 

Subroutine  SIZE  and  OBLATE  were  merged  together  to 
check  on  the  calculated  pressure  and  skin  drag  coefficient 
needed  for  the  six  equations  of  motion.  The  size  and  ter¬ 
minal  velocity  of  an  oblate  spheroid  is  fixed  only  after 
the  drag  coefficients  from  subroutine  OBLATE  are  numerically 
calculated.  From  Pitter  et  al . ^  >  35  ^  Table  o  .  2 . 1  shows  the 
comparison  between  results  reported  by  them  to  the  present 
research . 

Finally,  to  check  the  logic  of  NGCEFF  and  the  inte¬ 
gration  subroutine  COLL,  the  code  was  programmed  to  deter¬ 
mine  the  spherical  and  nonspherical  gravitational  collision 
efficiencies  for  oblate  spheroids  with  an  axis  ratio  equal 
to  0.999.  If  correctly  written,  the  results  of  efficiency 
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calculations  should  be  the  same  and  therefore  the  colli¬ 
sion  shape  factor,  8,  equal  to  one.  These  results  also 

afford  a  chance  to  compare  their  values  with  those  pub- 

1  s 

lished  by  Pertmer  and  Loyalka  .  For  the  NGCliFF  code 
operation,  the  containment  atmospheric  conditions  were 
taken  to  be: 

p^  =  1.29  x  10  ^  gm/cm^ 

-  4 

Pf  =  2.0  x  10  gm/cm-sec 

corresponding  to  a  temperature  of  120°F  and  pressure  of  0 
psig  for  air.  The  first  time  the  code  was  run  the  follow¬ 
ing  system  parameter  held  true:  AR  *  0.999,  Rej  *  0.49984, 
Cp1  =  52.261,  semi-major  axis  length  =  0.002574  cm,  p^  = 
2.27,  V°°  =  15.05  cm/sec,  v,  _  =  0.002574  cm.  The  values 

A  1  y  m 

of  the  small  particle  were:  r7  =  0.001  cm,  p7  =  2.27, 

Re2  =  0.0319,  V<»v  ^  ~  2.472  cm/sec.  System  parameters  were: 
Stk  =  2.044  ,  I.  =  2.277  ,  I  =  1.0,  y  =  1.0,  ii  =  1.0,  \<*>l  = 
16.38  cm/sec.  Results  of  NGCFFF  were: 

Volume  Basis  Mass  Basis 

eN  0.7102  0.7102 

eg  0.7099  0.7099 

6  =  1.000  (mass  ratio  <  0.059) 
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The  mass  ratio  is  important  to  check  since  values  greater 
than  0.34  will  tend  to  tilt  the  oblate  spheroid.  In  the 
present  case  it  is  not  important  since  the  larger  particle 
is  basically  a  sphere.  The  definitions  and  relationships 
for  the  parameters  listed  above  are  found  in  Chapters  IV  and 
V.  It  is  noted  that  in  general  the  terminal  velocity  of  the 
oblate  spheroid,  V°°  ,  is  not  equal  to  the  terminal  velocity 

a 

of  its  volume  equivalent  sphere,  V°°v  except  when  AR  = 

0.999.  Also,  V°°v  ^  is  not  equal  to  the  Stokes  terminal 
velocity,  V«=^,  based  on  the  oblate  spheroid's  volume  equiva¬ 
lent  radius,  because  the  former  ’ s  calculated  using  the  drag 
coefficient  from  the  Navier- Stokes  equation  whereas  the 
latter  makes  use  of  Stokes  terminal  velocity  expression 
(Re  =  0). 

For  a  similar  case  as  given  above  but  for  a  different 
small  aerosol  particle  with  a  =  0.0015  cm,  the  follow- 

ing  values  were  computed  by  NGCEFF :  Re2  =  0.107b,  V°°v  ^  = 
5.562,  Stk  =  10.35,  L  =  0.1153,  I  =  1.0,  y  =  1.0,  g  =  1.0, 

=  16.38  cm/sec, 

Volume  Basis  Mass  Basis 

eN  0.7866  0.7866 

es  0.7864  0.7864 

g  =  1.000  (mass  ratio  <  0.98). 
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A  third  case  was  run  for  an  oblate  spheroid  with  a 
semi-major  axis  length  equal  to  0.0038  cm.  The  following 
values  applied:  AR  =  0.999,  Re^  =  1.5,  =  2.27  g/cc, 

=  19.37,  V«a  =  30.2  cm/sec,  r^  m  =  0.003848  cm,  r^  y  = 
0.001,  p£  =  2.27  g/cc,  Re2  =  0.03189,  V°°2  =  2.47  cm/sec, 

Stk  =  1.106,  L  =  0.004566,  I  =  1.0,  y  =  1.0,  6  =  1.0,  V»1  = 
36.58  cm/sec.  The  results  were 

Volume  Basis  Mass  Basis 

eM  0.80787  0.80787 

N 

eg  0.80755  0.80755 

B  =  1.000  (mass  ratio  <  0.0176). 

Table  6.2.2  summarizes  these  results  and  compares  them 

1 2 

with  those  reported  by  Pertmer  .  Additional  cases  were 
available  for  large  particle  radius  38.5  pm  but  are  not 
shown  since  Pertmer's  results  for  particles  with  radii 
greater  than  30  pm  are  based  on  densities  less  than  2.27 
g/cc.  GCEFF  made  use  of  Figure  2. 1.2.1  for  establishing 
densities  of  colliding  particles. 

6.3  Results  from  the  NGCEFF  Code  -  B  Factors 


Collision  efficiencies  were  computed  for  oblate 
spheroids  of  axis  ratio  0.05  and  Reynolds  number  2.5  and 
spherical  particles  from  12.9  to  20.2  pm  in  radius.  Flow 
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Table  6.2.2 

Comparison  of  Calculated  Gravitational  Collision 
Efficiencies  for  Oblate  Spheroids  (Axis  Ratio  =  0.999) 

for  Several  Cases* 


Radii  (ym)  Ratio  Efficiency  Results 

r^y  r2>v  (r2>v/r1>y)  Present  Pertmerl* 


25.7 

10  0.39  0.71  0.66 

15  0.58  0.79  0.73 

38.5 

10  0.26  0.81  0.89** 


*From  Table  6. 3. 1.1,  Reference  12. 

**Extrapolated  Value. 


13b 


fields  around  the  oblate  spheroid  and  its  volume  equivalent 
sphere  were  taken  from  previous  computations  discussed 
earl ier . 

In  order  to  compute  the  collision  efficiencies  it  was 
necessary  to  determine  the  size  and  terminal  velocity  of 
the  oblate  spheroid.  Size  and  terminal  velocity  were  com¬ 
puted  j^sing  Equations  S.  1.1. 1.5  and  6. 1.1. 2.  Excellent 
agreement  was  found  between  the  terminal  velocities  calcu¬ 
lated  from  the  present  results  and  the  results  reported  by 
Pitter  et  al.^6. 

Table  6.3.1  and  6.3.2  lists  all  the  important  values 
and  quantities  for  the  case  being  considered.  Asterisks 
denote  input  values.  The  remaining  values  are  determined 
by  NGCEFF.  Additionally,  the  number  of  radial  and  angular 
steps  were  49  and  31,  respectively,  with  a  radial  increment 
size  of  0.1.  A  collision  was  scored  if  the  separation  dis¬ 
tance  was  less  than  10" 5  cm. 

Table  6.3.3  lists  the  collision  efficiencies  based  on 
geometric  sweepout,  i.e.. 


£g '  (r 


EP, 


+  r 


2,v 


(6.5.1) 


where  a  equals  semi-major  axis  length  and  the  remaining 
terms  are  defined  by  Equation  5.1.1.1.12.  These  were  cal¬ 
culated  for  comparison  purposes  with  those  reported  by 


Table  6.3.1 


Oblate  Spheroid  Values  and  System  Parameters  lor 
Collision  Results 


Oblate  Spheroid 

0.0116 

0.05* 

2.5* 

1.0* 

0.92* 

10,8 

19.4 

0.00426 

0.894 

31.03 

18.9 

0.00414 


Semi-Major  Axis  Length  (cm) 

Axis  Ratio 
Reynolds  Number 
Density  (g/cc) 

Density  Correction  Factor 
Drag  Coefficient 
Terminal  Velocity  (cm/sec) 

Volume  Equivalent  Radius  (cm) 
Equivalent  Reynolds  Number 
Equivalent  Drag  Coefficient 
Equivalent  Terminal  Velocity  (cm/sec) 
Mass  Equivalent  Radius  (cm) 


Collision  Parameters 


0.847 

0.00993 

1.00 

1.09 

1.00 

21.8 


Stokes  Number  (Stk) 

Interaction  Number  (L) 

Inertia  Number  (I) 

Ratio  of  Particle  Densities  (y) 

Ratio  of  Particle  Buoyancy  Effect  (8) 
Stokes  Terminal  Velocity  of  Volume 
Equivalent  Sphere 


Fluid  Parameters 


980.6* 

1.713  x  10 
6.500  x  10 
1.007  x  10 


-4* 
-  6  * 
-3* 


Gravitational  Constant  (cm/sec  ) 
Fluid  Viscosity  (gm/cm  x  sec) 
Mean  Free  Path  (A)  of  Fluid  (cm) 
Density  of  Fluid  (g/cc) 


DVOGER  Initial  Values 

1.0  x  10"!*  Initial  Time  Step  (sec) 
1.0  x  10\~*  Maximum  Time  Step  (sec) 
1.0  x  10"  Minimum  Time  Step  (sec) 


*Input  Values  to  NGCEFF 


Input  Values  to  NGCEFF 
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26 

Pitter  et  al .  .  This  was  not  possible,  however,  because 

the  annular  collision  domain  reported  by  Pitter  et  al .  was 
not  observed  in  any  of  the  present  cases.  As  discussed  in 
Section  6.1.2,  an  error  in  programming  may  have  resulted 
in  the  generation  of  the  annular  collision  domains  reported 
by  Pitter  et  al . . 

Table  6.3.3  lists  the  final  results.  These  are  signi¬ 
ficant  results  since  they  confirm  the  belief  that  it  should 
be  possible  to  theoretically  determine  nonspherical  colli¬ 
sion  efficiencies  and  the  collision  shape  factor  3.  By 
improving  on  the  methods  developed  in  this  research,  it  will 
be  possible  to  create  tables  of  shape  factors  for  various 
nonspherical  aerosols. 

Finally,  developers  of  aerosol  rate  codes  such  as 
HAARM-3  can  now  place  more  confidence  on  the  usage  of  a 
collision  shape  factors.  Current  practice  involves  running 
these  codes  with  different  3's  until  a  fairly  good  fit  to 
the  experimental  data  is  achieved.  By  using  shape  factors 
like  those  in  Table  6.3.4,  the  testing  of  these  codes  is 
more  rigorous  and  may  result  in  revealing  relationships 
that  were  being  masked  by  the  present  "curve  fitting" 
procedures . 


CHAPTER  VII 


CONCLUSIONS  AND  RECOMMENDATIONS 

7 . 1  Nonspherical  Gravitational  Collision  Efficiencies 

The  results  of  the  nonspherical  gravitational  colli¬ 
sion  efficiencies  indicate  that  it  is  inappropriate  to 
assume  that  agglomerates  collide  with  the  same  efficiencies 
as  spherical  particles.  During  the  early  stages  of  agglom¬ 
erate  growth,  there  will  be  collisions  between  chain-like 
particles  and  spherical  particles,  and  in  this  case  the 
results  from  Chapter  VI  show  clearly  that  the  collision 
shape  factor  8  is  necessary  to  modify  the  spherical  gravi¬ 
tational  collision  efficiencies.  For  the  case  considered 
in  Chapter  VI,  a  8  factor  equal  to  4.0  was  determined. 

The  abrupt  cutoffs  in  collision  efficiency  reported 

3  5 

by  Pitter  and  Pruppacher  was  not  observed  from  the  compu¬ 
tational  results.  As  discussed  in  Chapter  VI,  the  annular 
collision  domain  suggested  by  the  results  of  Pitter  and 
Pruppacher  is  most  likely  due  to  an  error  in  their  computer 
program.  Time  did  not  permit  a  change  in  their  program  to 
determine  if  different  results  would  be  calculated  but  when 
NCGEFF  was  modified  to  include  their  expressions  for  com¬ 
puting  the  collision  separation  distance,  some  of  their 
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results  were  obtained.  However,  considerable  problems 
occurred  because  the  separation  distances  were  always  being 
calculated  incorrectly,  which  resulted  in  the  spherical 
particle  inside  the  oblate  spheroid.  This  same  problem 
was  experienced  by  Pitter  and  Pruppacher.  When  the  correct 
expressions  were  substituted  into  NCGEFF,  this  problem  dis¬ 
appeared  . 

The  results  of  the  Navier-Stokes  solution  to  the  flow 
around  an  oblate  spheroid  showed  that  the  calculations  are 
very  sensitive  to  the  difference  expression  used  to  model 
the  boundary  conditions.  Higher  order  terms  [  >  0(h)] 
caused  the  solution  to  diverge;  the  smaller  the  error  assoc¬ 
iated  with  the  expression  was,  the  greater  the  divergence 
observed.  This  is  a  numerical  problem  with  the  finite  dif¬ 
ference  method,  although  it  was  not  suspected  until  after 
considerable  effort  had  been  spent  in  finding  the  part  of 
the  vorticity  field  which  caused  the  onset  of  divergence. 

The  use  of  cubic  splines  to  interpolate  the  velocity 
fields  associated  with  the  oblate  spheroid  was  superior 
when  compared  with  the  finite  difference  method.  The 
indeterminate  form  at  each  boundary  was  easily  calculated 
by  using  cubic  splines.  The  calculations  were  verified 
for  the  case  where  the  Reynolds  number  equals  zero  by  using 
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Oberbeck's  solution.  Boundary  calculations  were  further 
improved  by  using  Gregory-Newton  forward  and  backward  extra¬ 
polation  formula. 

7 . 2  Collision  Shape  Factor  (0) 

The  concept  of  the  collision  shape  factor  is  not  new 
but  in  this  dissertation  it  is  rigorously  derived  from  basic 
definitions.  Furthermore,  detail  information  in  Chapter  IV 
showed  that  the  gravitational  collision  efficiency  based  on 
mass  equivalent  radii  is  not,  in  general,  equal  to  the  gravi 
tational  collision  efficiency  based  on  volume  equivalent 
radii,  and  therefore  the  0  factor  will  be  incorrect  if  the 
same  basis  is  not  used  for  both  the  spherical  particle  and 
the  nonspherical  particle. 

The  calculation  of  6  factors  was  demonstrated  in 
Chapter  VI.  These  were  extremely  difficult  calculations  to 
do,  requiring  considerable  computer  time  and  hence  financial 
support  and  commitment  in  order  to  complete  this  research. 
However,  more  0  factors  are  necessary  if  any  confidence  is 
to  be  placed  in  their  usage  in  such  aerosol  rate  codes  like 
HAARM-3.  It  is  impossible  to  calculate  every  0  factor 

12 

needed  and  therefore  results  from  such  programs  like  GCEFF 
should  be  used  in  conjunction  with  a  set  of  0  factors. 
Currently  aerosol  rate  codes  make  use  of  the  concept  of  a 
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0  shape  factor  and  an  average  0  shape  factor  is  established, 
not  from  a  priori  knowledge,  but  from  curve  fitting  with 
experimental  data. 

7 . 3  Recommendations  for  Future  Work 

Basic  aerosol  research  is  applies  le  to  many  scientific 
fields  and  technologies.  The  research  reported  herein  is 
not  only  useful  to  the  developers  of  aerosol  rate  codes  for 
an  HCDA  in  an  LMFBR  but  equally  useful  to  developers  of  simi¬ 
lar  aerosol  rate  codes  for  postulated  generalized  working 
accidents  for  light  water  reactors  (LWR) ,  high  temperature 
gas-cooled  reactors  (HTGR)  and  gas-cooled  fast  reactors 
(GCFR) .  This  research  is  equally  important  to  researchers 
of  meteorologic  phenomena  and  environmental  pollution 
effects.  Thus,  basic  aerosol  research  involving  collisional 
mechanisms  impacts  many  areas  of  science  and  technology. 

In  this  section,  research  areas  where  future  work 
appears  worthwhile  are  outlined.  First  the  areas  in  which 
the  work  of  this  dissertation  can  be  improved  are  presented. 
Then  areas  not  considered  by  this  dissertation  that  need  to 
be  investigated  in  order  to  more  fully  understand  the  mecha¬ 
nisms  and  processes  that  are  influencing  particle  phenomena 
arc  highlighted. 


7.3.1  Supplementary  Dissertation  Research 


7. 3. 1.1  Tabulation  of  Collision  Shape  Factors 

It  was  noted  in  Chapter  VI  that  the  gravitational  col¬ 
lision  efficiencies  for  nonspherical  and  spherical  parti¬ 
cles  radically  differ  from  each  other.  Two  efforts  need 
undertaking.  First,  collision  shape  factors  for  oblate 
spheroids  and  spheres  for  various  densities  and  containment 
parameters  are  necessary  to  determine  bounds  on  this  very 
important  shape  factor.  Completion  of  this  phase  would  lead 
to  the  next  effort,  which  is  to  calculate  the  collision 
shape  factors  for  other  nonspiierical  particles  such  as  pro¬ 
late  spheroids.  When  completed  the  upper  and  lower  limits 
on  the  collision  shape  factor,  6,  would  be  better  defined 
(on  a  sound  theoretical,  basis)  then  it  had  been  before  in 
aerosol  science. 

7. 3. 1.2  Numerical  Methods 

Any  effort  to  tabulate  collision  shape  factors  needs 
to  address  improved  numerical  schemes  to  solve  the  Navier- 
Stokes  equation  and  the  six  linear  differential  equations 
of  aerosol  motion  and  interaction.  The  finite  element  method 
offers  promise  in  solving  the  Navier-Stokes  equation  in  both 
spherical  and  spheroids  coordinate  systems.  Likewise,  addi¬ 
tional  research  of  the  stiff  behavior  of  the  six  equations 
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of  aerosol  interaction  may  produce  a  more  efficient  inte¬ 
gration  routine.  These  areas  of  research,  if  successful, 
will  have  a  significant  impact  on  reducing  computational 
expenses  associated  with  shape  factor  tabulation. 

7. 3. 1.3  NGCEFF  Code  Modification 

The  computer  code  NGCEFF  listed  in  Appendix  1  to  deter¬ 
mine  the  collision  shape  factor,  3,  is  a  very  general  code. 

In  its  present  configuration  it  can  determine  aerosol  parti¬ 
cle  trajectories,  drag  forces  between  particles  and  velocity 
fields  around  oblate  spheroids.  This  information  is  of 
scientific  interest  and  should  be  merged  with  graphics  pack¬ 
ages  for  display  and  tabulated  for  future  reference.  NGCEFF 
can  be  modified  or  expanded  as  required  without  a  major 
rewriting  to  incorporate  new  information  such  as  kinetic 
corrections  to  drag  forces. 

7. 3.1.4  Functional  Representation  of  the  Collision  Shape 
Factor 

Since  the  collision  shape  factor,  3,  plays  an  impor¬ 
tant  role  in  aerosol  rate  code  such  as  HAARM-3,  effort  should 
go  into  the  development  of  an  efficient  and  accurate  func¬ 
tional  representation  of  3  after  tabulation  of  its  values  is 
accomplished.  In  the  past  the  data  were  fit  using  a  natural 


cubic  spline  to  interpolate  between  known  values.  This  may 
not  be  the  best  way  to  represent  6,  because  of  the  increase 
in  the  number  of  variables  and  therefore  research  is  nec¬ 
essary  to  define  the  problem  and  formulate  an  approximate 
solution. 

7.3.2  Additional  Research  Topics 

7. 3. 2.1  Synergistic  Effects 

It  is  stipulated  that  the  three  coagulation  mechanisms 
are  separative  and  additive  for  the  collision  kernel,  <p , 
Equation  2.3.1.  Possible  synergistic  effects  of  one  mecha¬ 
nism  on  another  should  be  investigated  in  order  to  verify 
this  assumption  or  determine  what  corrections  are  necessary 
to  model  observed  phenomena  and  its  associated  theory. 

7. 3. 2. 2  Coalescent  Effects 

One  effect  that  needs  thorough  investigation  is 
coalescence  of  particles.  This  effect  must  be  studied  both 
by  the  experimentist  and  theorist  and  a  method  for  calcu¬ 
lating  the  fraction  of  particles  that  coalesce  determined. 

7. 3. 2. 3  Knudsen  Drag  Forces 

As  reported  in  Appendix  3,  Knudsen  drag  forces  can  sig¬ 
nificantly  improve  the  accuracy  of  the  computational  method 
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used  to  determine  the  collision  efficiency.  Considerable 
more  effort  is  needed  in  this  area  even  if  it  is  extremely 
complicated.  This  research,  if  successful,  would  be  extremely 
important  to  the  developers  of  aerosol  rate  codes  and  to  the 
scientific  community  as  a  whole. 

7. 3. 2.4  Other  Aerosol  Effects 

Many  other  aerosol  effects  must  be  taken  into  consider¬ 
ation.  Some  of  these  include  radiation  and  electrical 
charge  effects,  water  vapor  (humidity)  and  thermal  gradient 
effects,  and  dif fusiophoretic ,  inertial  and  turbulent  deposi¬ 
tion  effects. 

The  complete  description  of  the  processes  affecting  an 
aerosol  is  extremely  complicated.  It  is  important,  however, 
to  concentrate  on  understanding  these  processes  individually 
and  then  together  before  any  real  confidence  can  be  placed 
on  codes  like  HAARM-3. 

7. 3. 2. 5  Experimental  Data 

It  is  important  to  compare  the  results  of  the  present 
research  with  experimental  data.  An  important  consideration 
here  is  to  study  the  effects  of  particle  shape  on  the  gravi¬ 
tational  collision  efficiency  and  collision  shape  factor. 
Aerosol  characteristics  and  fluid  parameter  of  post-HCDA  con¬ 
tainments  need  to  be  included  into  the  design  of  such  experi¬ 


ments  . 
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APPENDIX  1 


COMPUTER  CODE  LISTINGS 

The  listings  of  the  computer  code  NCGEFF  is  contained 
in  this  Appendix.  The  complete  descriptions  of  the  code  and 
of  code  usage  can  be  found  in  Chapter  V. 

The  code  utilizes  the  following  subroutines  from  the 
International  Mathematical  and  Statistical  Libraries,  Inc. 
(IMSL) :  ZFALSE ,  ICSICU,  and  ICSEVU.  The  IMSL  subroutines 

are  not  listed.  Calculation  of  the  eigenvalues  of  the  two 
matrices  associated  with  the  vorticity  transport  equations 
is  done  with  the  UNPACK  routines  and  are  available  from 
Argonne  Laboratory  or  IMSL.  Subroutine  DVOGER  has  been  merged 
with  NCGEFF  since  it  is  no  longer  a  current  IMSL  subroutine. 

It  was  utilized  in  lieu  of  DGEAR  because  of  certain  features 
not  available  in  DGEAR. 
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APPENDIX  2 


THREE  FUNCTION  ROUTINES  FOR  SPHERICAL  PARTICLES 

Computer  listings  of  three  FORTRAN  IV  function  subpro¬ 
grams  are  included  in  this  appendix.  The  spherical  gravita¬ 
tional  collision  efficiency.  Eg,  for  particles  with  constant 

3  3  3 

densities  of  0.30  g/cm  ,  1.00  g/cm  or  2.27  g/cm  can  be 
generated  by  using  these  subprograms  in  conjunction  with 
two  IMSL  double  precision  routines,  ICSICU  and  ICSEVU.  Me 
A2.1  summarizes  the  key  parameter  associated  with  these 
tines. 

All  values  of  Eg  were  determined  by  means  of  the  GCEFF 

12  12 
code  and  associated  programs  called  SPLFIT  and  GPROC. 

GPROC  was  a  temporarily  created  TSO  executive  program  written 

by  the  author  to  handle  the  data  from  GCEFF  and  SPLFIT  and 

to  create  the  necessary  DATA  statement  for  the  GEPS  routines. 

Since  it  made  use  of  UMC  installation  depended  commands,  it 

was  not  included  in  this  dissertation.  The  author  wishes  to 

express  his  appreciation  to  Mr.  Takoshi  Enomoto  for  his 

assistance  in  running  GCEFF,  SPLFIT,  and  GPROC. 
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Table  A2.1 

GEPS  FUNCTION  CAPABILITIES  AND  LIMITATIONS 


— 

CALL  NAME (Arguments) 

Density 

g/cm3 

Large  Particle 
Radius , cm* 

Small  Particle 
Radius , +cm 

GEPS1 (A1 , A2) 

2.27 

10~4<A1<10~2 

A2 

GEPS2 (A1 , A2) 

1.00 

10‘4<a1<10'2 

A2 

GEPS3 (A1 , A2) 

0.30 

10'4<Al<10'2 

A2 

*When  the  ratio  of  the  smaller  particle  to  the  larger 
particle  is  less  than  0.1  or  greater  than  0.9,  or  the  large 

A  1 

particle  radius  is  less  than  10"  cm  or  greater  than  10  cm, 
Fuchs  expression  is  used  by  the  GEPS  routines.  Fuchs  expres¬ 
sion  is: 


eS  =  1,50  (A1  ♦  A*) 

A2  is  determined  by  (A)(A1)  where  A  is  A2/A1  and  is 
bound  as : 


0.1  <  A  <  0.9. 
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APPENDIX  3 


KINETIC  CORRECTIONS  TO  AEROSOL 
GRAVITATIONAL  COLLISIONAL  EFFICIENCY 

Sensitivity  analyses  have  shown  that  the  gravitational 
collision  efficiency  influences  post-hypothetical  core  dis¬ 
ruptive  accidents  aerosol  behavior  in  LMFBR  containment  in 

17  18 

important  ways  .  It  was  noted  by  Pertmer  and  Loyalka 

that  the  assumptions  used  in  computing  two  body  drag  forces 
influences  the  calculated  collisional  efficiency  strongly. 

In  fact,  the  Stokes  and  Oseen  approximations  to  drag  forces 
can  lead  to  results  for  collisional  efficiency  that  are  dif¬ 
ferent  from  each  other  by  orders  of  magnitude.  In  their 
work  it  was  also  noted  that  the  use  of  Oseen  drag  forces  pro 
vided  results  that  agreed  well  with  some  available  experi¬ 
mental  data  for  large  particles  (radius  >  30  pm) .  For  small 
particles  (radius  <  30  ym)  one  would  generally  expect  Stokes 
to  provide  better  agreement  with  experimental  data,  but 
Pertmer  and  Loyalka  found  that  Stokes'  results  differed 
from  experimental  data  by  about  two  orders  of  magnitude. 
Curiously,  the  Oseen  results  showed  less  deviations  from 
experimental  data  (at  most,  one  order  of  magnitude;  see 
Figure  A3. 1) . 
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This  circumstance  led  to  a  reexamination  of  the  Stokes' 
results.  The  inadequacy  of  the  Stokes’  drag  forces  for 
small  particles  and  small  separations  between  otherwise 
large  particles  was  discussed  in  Chapter  III.  An  acceptable 
analysis  of  the  forces  between  small  particles  (<  1  y)  or 
large  Knudsen  numbers  (kn  >0.1)  requires  some  type  of  ki¬ 
netic  corrections,  i.e.,  requires  that  the  Knudsen  drag 
forces  on  the  two  spherical  particles  be  obtained  by  solving 
the  relevant  boundary  value  problem  in  the  framework  of  the 
kinetic  theory  of  gases.  Pertmer  and  Loyalka  worked  out 
the  mathematics  necessary  for  the  general  case  of  two  spher¬ 
ical  particles  and  the  details  can  be  found  in  References  12 
and  46. 

Basically,  the  Knudsen  drag  forces  are  determined  by 
solving  the  relevant  boundary  value  problem  for  particles 
interacting  in  a  monatomic  simple  gas.  This  is  a  difficult 
problem  and  research  is  currently  underway  to  obtain  expres¬ 
sions  for  such  forces.  Meanwhile,  useful  results  can  be  ob¬ 
tained  by  using  recently  determined^  velocity  profiles 
for  a  sphere  moving  in  a  rarefied  gas  and  the  method  of  super¬ 
position.  The  work  of  Lee  and  Loyalka  provides  numerical 
results  for  velocity  profiles  for  all  kn  by  the  use  of  the 
BGK  model.  Tomoeda  also  considered  the  BGK  model  and  ob¬ 
tained  results  for  kn  <  .1  by  using  a  technique  due  to  Sone. 
For  kn  <  .1,  the  results  of  Lee  and  Loyalka  are  in  agreement 
with  those  of  Tomoeda.  Using  the  superposition  principle, 
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Knudsen  and  Stokes  drag  forces  can  now  be  approximated  so 
that  the  six  nondimens ional ,  linear  first-order  differential 
equations  of  motion  can  be  solved  to  calculate  the  collisional 
efficiencies  for  spherical  particles.  Pertmer's  code, 

GCEFF,  incorporated  all  necessary  features  except  for  the 
Knudsen  velocity  profiles.  It  was  therefore  decided  to  de¬ 
rive  an  approximation  to  Knudsen  velocity  profiles  which  were 
asymptotically  matched  to  Stokes  at  a  few  mean  free  paths 
from  a  sphere,  and  to  incorporate  the  results  into  GCEFF. 

A  convenient  representation  of  the  velocity  profiles  for 
kn  <<  1  can  be  obtained  by  using  Tomoeda's  results.  If  we 
consider  spheres  of  sufficiently  large  size  such  that  kn  <<  1, 
then  to  a  good  degree  of  approximation,  we  can  write: 


U-  — 
.  xjt 
COSlj) 


■!<£>«-  i  <T->  -  JT 


[1.0161  - 


1  {_L)  .  (!i,2 

2  la.;  1  r  ' 


(1.0161  - 


0.072488 


)]> 


+  fCT^T^Tj)  (A3.1) 

-  -  J  <T->  tl-°  -  C1-0161  -  7  £ 

+  (^  -  (1.0161  -  0.32644  [i-])  i- 

•  i~)Z)  *  f'(T0 


(A3. 2) 
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fdj.Tj.Tj)  =  3.0(i-)2  {-0.653551  Tj  (-^3 

r"ai 

+  1.471319  T2(— 


-  1.117159  T3(— j-i) 

f' (T0,T1,T2,T3)  =  3  (-0.653551  TQ(— j~) 

r-  a- 

+  1.471319  T^— x^) 
r-a- 

-  1.117159  T2(— j—)) 

3  A  2  r~a- 

-  f  (~)Z  {>1.01618  TjC-jpi) 

r'ai 

+  4.517256  T2(-ri) 


(A3. 3) 


where 


r>  a. 

-  2.504171  T3(-~)> 


Val 


■  tan **y/x 


and  these  quantities  are  defined  as 


(A3. 4) 


(A3. 5) 
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^  aerosol  particle  of  radius  a^,  i  =  1,2 

r  radial  distance  from  center  of  particle 
X  a,  mean  free  path  for  containment  gas 
y  ^  Cartesian  coordinate  distance  from  center  of 
particle,  y  *  ?  =>  sin$  =  ±  1 
x  Cartesian  coordinate  distance  from  center  of 
particle,  x  *  ?  *>  sin$  =  0 
TM(arg)  'x*  Abramowitz  Functions46 


In  Table  A3.1  the  calculated  collisional  efficiencies 
for  larger  particles  radius  •  0.1  to  30  yin,  and  the  par¬ 
ticle  radii  ratio  a  *  *2^1  “  to  0.90  is  compared  with 

Stokes'  superposition  results  from  Pertmer.  Considerable 
discrepancy  is  noted  for  ratios  <  0.40  and  for  collisions 
between  all  size  particles  <  5  ym.  Note  that  Stokes'  results 
always  lead  to  collisional  efficiencies  that  are  too  low, 
and  hence  use  (in  aerosol  behavior  codes)  of  collisional 
efficiencies  based  on  Stokes  drag  forces  will  lead  to  very 
conservative  predictions  of  radioactive  release. 

Referencing  the  work  of  Tu  and  Shaw,48  the  modified 
GCF.FF  code  was  programmed  to  calculate  the  collisional  ef¬ 
ficiencies  between  drops  of  E.  Coli  and  water,  the  same  ex¬ 
perimental  parameter  used  in  the  work  of  Tu  and  Shaw.  Thus 

=  1.004  g/cm3 
P2  -  1.0  g/cm3 


*Upper  Value:  Stokes  with  Kinetic  Corrections 
Lower  Value:  Stokes  Only 


2SS 


pf  -  1.007  x  IQ*3  g/cm3 
X  -  7.37  x  10'6  cm3 
Vjc  “  1.713  x  10**  poise 

where  is  the  density  of  E.  Coli  droplets  and  P2  is  the 
density  of  water  droplets,  p^  is  the  density  of  air,  X  is 
the  mean  free  path  of  air,  and  y  is  its  viscosity.  The  pre¬ 
sent  results  together  with  the  experimental  data  and  several 
other  theoretical  results  are  given  in  Figure  A3.1.  It  is 
noted  that  the  present  results  describe  the  experimental 
data  well  for  >0.4,  but  for  i  0.40  the  present  results  do 
not  agree  with  experimental  data.  There  is  a  dearth  of  clear 
experimental  data  in  this  regime  and  therefore  it  is  not 
clear  if  the  disagreement  is  due  to  inadequacies  of  the  super¬ 
position  method  or  the  experimental  data  or  both.  The  pre¬ 
sent  work,  however,  emphasizes  the  improvement  offered  by 
modeling  the  effects  of  slip  and  this  justifies  further  ex¬ 
perimental  work  and  still  more  improved  theoretical  calcu¬ 
lations. 
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